Penyelesaian Persamaan Diferensial Ordiner dengan MATLAB Pada pemodelan kasus-kasus di teknik kimia seringkali dijumpai hasil persamaan matematikanya matematikanya berupa persamaan diferensial ordiner (PDO). Dalam Persamaan diferensial ordiner ini hanya dijumpai satu variabel bebas. Terdapat dua jenis PDO yang sering dijumpai dalam pemodelan matematis, yaitu: 1. PDO jenis initial value problem Kondisi batas yang diketahui nilainya hanya pada satu titik. 2. PDO jenis boundary value problem Kondisi batas yang diketahui nilainya pada lebih dari satu titik.
1. Persamaan Diferensial Ordiner Jenis Initial Value Problem (IVP) 1.1.
Penyelesaian dengan integrasi numeris
Dalam beberapa kasus pemodelan, PDO yang diperoleh bisa langsung diselesaikan dengan integrasi. Integrasi ini bisa dilakukan secara analitis maupun secara numeris. Jika integrasi dilakukan secara numeris, dalam MATLAB bisa dipergunakan fungsi INTEGRAL. INTEGRAL. Contoh 1. Dari pemodelan kasus 2 (Absorpsi gas dalam solven secara batch) diperoleh PDO sebagai berikut:
( = ∙ )
(co.1.1)
Dengan kondisi batas: t = 0; y = 0. Persamaan (1) di atas bisa diintegrasikan untuk mendapatkan mendapatkan hubungan y sebagai fungsi t sebagai berikut:
∫ = ∙ ∫ (−) ∙ = ∙ ∫ (−)
(co.1.2)
Waktu (t) yang diperlukan sehingga kadar gas yang keluar menjadi sebesar yout dapat dihitung dengan fungsi INTEGRAL dalam MATLAB. Kode pemrograman dalam MATLAB function contoh01PEMAT function contoh01PEMAT % Kasus 02 Absorpsi gas dalam solven secara batch % Penyelesaian dengan integrasi numeris clc clear
% Data untuk perhitungan polutan = 'CO2' 'CO2'; ; m = 1000; % jumlah solven dalam tangki, mol G = 1; % laju alir gas yang digelembungkan dalam tangki, mol/detik H = 0.186; % Konstanta Henry gas CO2 yin = 0.5; % fraksi mol CO2 dalam gas masuk % Menghitung waktu yang diperlukan supaya kadar gas keluar sistem menjadi % 10% kadar masuknya yout = 0.1*yin; % Kadar CO2 keluar sistem % Perhitungan secara numeris dengan fungsi INTEGRAL tnum = m*H/G*integral(@fungsi,0,yout); %Perhitungan secara analitis texact = m*H/G*log(yin/(yin-yout)); % Menampilkan hasil perhitungan fprintf('Kasus fprintf('Kasus 02 Absorpsi gas dalam solven secara batch.\n') batch.\n' ) fprintf('------------------------------------------------\n' fprintf('------------------------------------------------\n' ) fprintf('Gas fprintf('Gas %s berkadar %4.2f diserap dalam air.\n',polutan,yin) air.\n',polutan,yin) fprintf('Kadar fprintf('Kadar %s maksimum yang keluar adalah %4.2f.\n',polutan,yout) %4.2f.\n',polutan,yout) fprintf('Waktu fprintf('Waktu yang diperlukan:\n') diperlukan:\n') fprintf('\tPerhitungan fprintf('\tPerhitungan secara numeris %4.2f detik. \n',tnum) \n',tnum) fprintf('\tPerhitungan fprintf('\tPerhitungan secara analitis %4.2f detik.\n',texact) detik.\n',texact) function f = fungsi(y) function f f = 1./(yin-y); end end
Hasil perhitungan: Kasus 02 Absorpsi gas dalam solven secara batch. -----------------------------------------------Gas CO2 berkadar 0.50 diserap dalam air. Kadar CO2 maksimum yang keluar adalah 0.05. Waktu yang diperlukan: Perhitungan secara numeris 19.60 detik. Perhitungan secara analitis 19.60 detik.
1.2.
Penyelesaian dengan ODE solver
Jika hasil pemodelan matematis yang diperoleh berupa PDO, maka dapat diselesaikan dengan fungsi ODE solver yang ada di MATLAB. Ada beberapa jenis ODE solver dalam MATLAB, yaitu: Tabel 1. Berbagai fungsi ODE solver dalam MATLAB MATLAB Solver
Solves These Kinds of Problems
Method
ode45
Nonstiff differential equations
Runge-Kutta
ode23
Nonstiff differential equations
Runge-Kutta
Solver
Solves These Kinds of Problems
Method
ode113
Nonstiff differential equations
Adams
ode15s
Stiff differential equations and DAEs
NDFs (BDFs)
ode23s
Stiff differential equations
Rosenbrock
ode23t
Moderately stiff differential equations and DAEs DAEs Trapezoidal rule rule
ode23tb Stiff differential equations
TR-BDF2
ode15i
BDFs
Fully implicit differential equations
Semuanya mempunyai syntax yang sama dalam pemakaiannya. Fungsi yang biasanya dipergunakan adalah ODE45. Jika fungsi ODE45 terlalu lambat dalam menyelesaikan PDO atau diketahui PDO mempunyai hasil penyelesaian yang berubah tajam (stiff ( stiff ) maka lebih baik dipergunakan fungsi ODE15S. Contoh 2. Pemodelan Kasus 1: Pencampuran dalam tangki. PDO yang diperoleh:
= ( = )
(co.2.1) (co.2.2)
Masa larutan dalam tangki dan kadar garamnya setiap saat dapat dihitung dengan fungsi ODE45 dalam MATLAB. Kode pemrograman dalam MATLAB function contoh02aPEMAT function contoh02aPEMAT %Kasus 1: Pencampuran dalam tangki %Penyelesaian dengan ODE45 clc; clear % Data V0 = 100; x0 = 0; % F = 15; % L = 10; % xF = 1; %
% Massa air dalam tangki semula (kg) Kadar garam dalam tangki semula (kg/L) Debit larutan garam masuk ke tangki (kg/men) Debit larutan garam keluar dari tangki (kg/men) Kadar garam dalam larutan masuk ke tangki (kg/L)
% Penyelesaian dengan fungsi ODE45 tf = 10; % Waktu final (men) IC = [V0 x0]; [t Vx] = ode45(@ode_kasus01,[0 tf], IC); V = Vx(:,1); % Massa larutan dalam tangki x = Vx(:,2); % Kadar garam larutan % Menampilkan hasil perhitungan dalam bentuk grafik
subplot(1,2,1) plot(t,V) title('V title('V vs t') t') xlabel('Waktu, xlabel('Waktu, men'); men'); ylabel('Massa ylabel('Massa larutan dalam tangki, kg') kg') subplot(1,2,2) plot(t,x) title('x title('x vs t') t') xlabel('Waktu, xlabel('Waktu, men'); men'); ylabel('kadar ylabel('kadar garam larutan dalam tangki, kg') kg') function dVxdt = ode_kasus01(tfun, Vxfun) function dVxdt dVxdt = zeros(2,1); Vfun = Vxfun(1); xfun = Vxfun(2); dVxdt(1) = F - L; dVxdt(2) = F/Vfun * (xF - xfun); end end
Hasil perhitungan:
Contoh 3. Perbaikan pemrograman pemrograman komputer untuk contoh 2. Ada beberapa kondisi dalam contoh 2 yang harus diantisipasi dalam program komputer yang dibuat supaya hasil perhitungan yang diperoleh sesuai dengan kondisi sebenarnya. Beberapa kondisi tersebut adalah:
= 0 dan = 0 . Jika L < F, suatu saat tangki akan meluap setelah V > Vmax. J ika V ≥ Vmax maka = 0 dan tetap sesuai persamaan co.2.2.
1. Jika L > F, suatu saat tangki akan kering (V = 0). Jika J ika V = 0 maka 2.
Dengan kata lain, persamaan co.2.1. hanya berlaku jika 0 ≤ V ≤ Vmax dan co.2.2. hanya berlaku jika V ≥ 0. Kode permrogramannya dalam MATLAB: function contoh03aPEMAT function contoh03aPEMAT %Kasus 1: Pencampuran dalam tangki %Penyelesaian dengan ODE45 clc; clear % Data V0 = 100; % Massa air dalam tangki semula (kg) x0 = 0; % Kadar garam dalam tangki semula (kg/L) F = 15; % Debit larutan garam masuk ke tangki (kg/men) L = 10; % Debit larutan garam keluar dari tangki (kg/men) xF = 1; % Kadar garam dalam larutan masuk ke tangki (kg/L) Vmax = 150; % Massa maksimum yang dapat ditampung tangki % Penyelesaian dengan fungsi ODE45 tf = 50; % Waktu final (men) IC = [V0 x0]; [t Vx] = ode45(@ode_kasus01,[0 tf], IC); V = Vx(:,1); % Massa larutan dalam tangki x = Vx(:,2); % Kadar garam larutan % Menampilkan hasil perhitungan dalam bentuk grafik subplot(1,2,1) plot(t,V) title('V title('V vs t') t') xlabel('Waktu, xlabel('Waktu, men'); men'); ylabel('Massa ylabel('Massa larutan dalam tangki, kg') kg') subplot(1,2,2) plot(t,x) title('x title('x vs t') t') xlabel('Waktu, xlabel('Waktu, men'); men'); ylabel('kadar ylabel('kadar garam larutan dalam tangki, kg') kg') function dVxdt = ode_kasus01(tfun, Vxfun) function dVxdt dVxdt = zeros(2,1); Vfun = Vxfun(1); xfun = Vxfun(2); if Vfun if Vfun >= 0 dVxdt(1) = F - L; dVxdt(2) = F/Vfun * (xF - xfun); end if Vfun if Vfun >= Vmax dVxdt(1) = 0; end end end
Hasil perhitungan:
Dalam grafik di atas untuk V vs t ada masalah yaitu: Di sini grafik V melebihi nilai maksimum yang diijinkan (Vmax). Juga terlihat grafiknya tidak mendatar namun ada sedikit fluktuasi. Hal tersebut tidak sesuai dengan kondisi sebenarnya.
Kesalahan di atas dapat terjadi karena langkah perhitungan maksimum (maximum step size; MAXSTEP) yang dipergunakan dalam ODE45 terlalu besar. Perlu dilakukan modifikasi dalam perintah ODE45 sebagai berikut: % Penyelesaian dengan fungsi ODE45 tf = 100; % Waktu final (men) IC = [V0 x0]; pilihan = odeset('MaxStep' odeset('MaxStep',0.01); ,0.01); [t Vx] = ode45(@ode_kasus01,[0 tf], IC, pilihan);
Kode pemrograman yang lain sama dengan dalam file contoh03aPEMAT perhitungan menjadi:
di atas. Hasil
Keterangan lebih detail dapat dibaca dalam HELP pada MATLAB. 2. Persamaan Diferensial Ordiner Jenis Boundary Value Problem (BVP) PDO jenis BVP ini mempunyai lebih dari satu kondisi batas yang semuanya harus terpenuhi. Dalam MATLAB dapat dipergunakan fungsi BVP4C atau BVP5C untuk penyelesaian PDO jenis BVP. BVP4C dan BVP5C menyelesaikan PDO jenis BVP dalam bentuk:
0 ( ) = , dengan kondisi batas [ ] = [ ] 0
Contoh 4. Kasus 6. Distribusi Suhu pada Batang di Udara PD yang diperoleh:
4∙ℎ ( = ∙ )
Kondisi batas: a. Pada x = 0; T = Ts b. Pada x = L ada 3 kemungkinan: 1. T = Tf (konstan) 2. 3.
= 0 ℎ( = )
(co.4.1.) (co.4.2.) (co.4.3.1.) (co.4.3.2.) (co.4.3.3.)
Penyelesaian: PDO order 2 pada persamaan co.4.1. diubah menjadi PDO order 1 dengan cara sebagai berikut:
= 2 4∙ℎ ( ) = ∙ 2 )
Kondisi batas menjadi: a. Pada x = a; y2(a) = Ts b. Pada x = b ada 3 kemungkinan: kemungkinan: 1. y1(b) = Tf (konstan) 2. 2 3.
() = 0 2() = ℎ (1 )
(co.4.4) (co.4.5) (co.4.6.) (co.4.7.1.) (co.4.7.2.) (co.4.7.3.)
Contoh 4.1. Kode MATLAB jika kedua ujung batang bersuhu konstan (persamaan co.4.7.1.). function contoh4aPEMAT function contoh4aPEMAT % Kasus 6. Distribusi Suhu pada Batang di Udara % Contoh penyelesaian PDO BVP clc; clear % Data untuk perhitungan k = 0.5; %konduktivitas panas (cal/det/cm/C) h = 0.0007; %koefisien perpindahan panas batang-udara (cal/det/cm2/C) L = 100; %panjang batang (cm) d = 1; %diameter batang (cm) Ts = 400; %Suhu permukaan luar pipa (C) Tu = 30; %Suhu udara (C) Tf = 100; % Suhu di ujung luar batang % Menghitung distribusi suhu dalam batang dengan fungsi BVP4C solinit = bvpinit(linspace(0,L,100),[Ts 0]); % Nilai tebakan y1 dan y2 sol = bvp4c(@odefun,@bcfun,solinit); % Penyelesaian untuk y1 dan y2 % x y T
Menampilkan nilai distribusi suhu (T) = linspace(0,L); = deval(sol,x); = y(1,:);
% membuat grafik hasil perhitungan plot(x,T) title({'Distribusi title({'Distribusi Suhu Dalam Batang'; Batang';... 'jika kedua ujungnya bersuhu konstan'}) konstan'}) xlabel('Panjang xlabel('Panjang batang, cm') cm') ylabel('Suhu ylabel('Suhu \circC') \circC') function dydx = odefun(x,y) function dydx dydx = zeros(2,1); dydx(1) = y(2); dydx(2) = 4*h/k/d*(y(1) - Tu); end function bc function bc = bcfun(ya,yb)
bc = zeros(2,1); bc(1) = ya(1)-Ts; bc(2) = yb(1)-Tf; end end
Hasil perhitungan:
Contoh 4.2. 4.2. Kode MATLAB jika salah satu ujung batang diisolasi (persamaan co.4.7.2.). co.4.7.2.). Kodenya sama dengan kasus 1, hanya fungsi untuk kondisi batasnya diubah menjadi: function bc = bcfun(ya,yb) function bc bc = zeros(2,1); bc(1) = ya(1)-Ts; bc(2) = yb(2); end
Hasil perhitungannya: perhitungannya:
Contoh 4.3. Kode MATLAB jika salah satu ujung batang dibiarkan tanpa isolasi (persamaan co.4.7.3.). Kodenya sama dengan kasus 1, hanya fungsi untuk kondisi batasnya diubah menjadi: function bc = bcfun(ya,yb) function bc bc = zeros(2,1); bc(1) = ya(1)-Ts; bc(2) = yb(2)+h/k*(yb(1)-Tu); end
Hasil perhitungannya: perhitungannya:
