Nilai Awal (IVP = initial value problems)
1. Persamaan diferensial biasa (PDB) orde satu tunggal
(Single, first order ordinary differential equation (ODE))
2. Sistem persamaan diferensial biasa
3. Persamaan diferensial parsial (PDE)
I. PDB orde satu tunggal
[Link] Euler Eksplisit
[Link] Euler Implisit Dapat diselesaikan dengan deret Taylor
[Link] Trapesium
[Link] Runge-Kutta Orde-4
1
Analitik :
dy dy
= x 2 . y ========> = x 2 dx
dx y
1 3
ln y = . x + C
3
x = 0 =====> y = 1 ====> C = 0
1 3
ln y = .x
3
1 x3
y = e3
x = 1 ======> y = 1,3956
dy
= x2 + y ======> Sulit untuk menyelesaikan secara analitik
dx
Solusi : penyelesaian secara numerik
2
Metoda Euler Eksplisit
Rumus umum :
yi+1 = yi + Δx.f(xi,yi) (i dimulai dari 0)
atau : yi = yi-1 + Δx.f(xi-1,yi-1) (i dimulai dari 1)
Contoh :
Tentukan nilai y pada x = 1 dengan menggunakan metoda euler eksplisit pada persamaan
berikut :
dy
= x2 y ; dimana y = 1 pada x = 0
dx
Penyelesaian :
Metoda euler eksplisit : yi = yi-1 + Δx * xi-1 * xi-1*yi-1
Pilih : Δx = 0,1 untuk xo = 0 dan yo = 1
i x yi = yi-1 + Δx * xi-1 * xi-1*yi-1 yi
1 0,1 y1 = yo + Δx * xo* xo*yo = 1 + 0,1. 02.1 1
2 0,2 y2 = y1 + Δx * x1* x1*y1 = 1 + 0,1. 0,12.1 1,001
3 0,3 y3 = y2 + Δx * x2* x2*y2 = 1,001 + 0,1. 0,22.1,001 1,005
4 0,4 y4 = y3 + Δx * x3* x3*y3 = 1,005 + 0,1. 0,32.1,005 1,014
5 0,5 y5 1,0303
6 0,6 y6 1,0560
7 0,7 y7 1,0940
8 0,8 y8 1,1477
9 0,9 y9 1,2211
10 1,0 y10 1,3200
Bandingkan jika :
Δx 0,05 0,02 0,01 0,0001 0,000001 analitik
x 1,35588 1,37922 1,38733 1,39553 1,39561 1,39561
3
clc
% Program : euler.m
h = 0.1; % h = step/interval
x = [0 : h : 1];
y = 0*x;
y(1) = 1;
for i = 2 : max(size(x))
y(i) = y(i-1) + h*x(i-1)*x(i-1)*y(i-1);
end
analitik = exp(1/3*x.^3);
disp (' ')
disp ('Hasil perhitungan :')
disp ('---------------------------------')
disp (' x y analitik')
disp ('---------------------------------')
for i = 1:max(size(x))
fprintf('%6.1f%10.5f%12.5f\n',x(i),y(i),analitik(i))
end
disp ('---------------------------------')
» euler
Hasil perhitungan :
---------------------------------
x y analitik
---------------------------------
0.0 1.00000 1.00000
0.1 1.00000 1.00033
0.2 1.00100 1.00267
0.3 1.00500 1.00904
0.4 1.01405 1.02156
0.5 1.03027 1.04255
0.6 1.05603 1.07466
0.7 1.09405 1.12113
0.8 1.14766 1.18610
0.9 1.22111 1.27507
1.0 1.32002 1.39561
---------------------------------
4
function F = f803(x,y)
F = x*x*y;
» [x y] = ode23('f803',[0:0.1:1],1);
» [x y]
ans =
0 1.0000
0.1000 1.0003
0.2000 1.0027
0.3000 1.0090
0.4000 1.0216
0.5000 1.0425
0.6000 1.0747
0.7000 1.1211
0.8000 1.1861
0.9000 1.2751
1.0000 1.3956
5
clc
% Program : rk.m
h = 0.1; % h = step/interval
x = [0 : h : 1];
y = 0*x;
y(1) = 1;
for i = 2 : max(size(x))
k1 = (x(i-1)^2)*y(i-1);
k2 = ((x(i-1)+h/2)^2)*(y(i-1)+h*k1/2);
k3 = ((x(i-1)+h/2)^2)*(y(i-1)+h*k2/2);
k4 = ((x(i-1)+h)^2)*(y(i-1)+h*k3);
y(i) = y(i-1) + (h/6)*(k1+2*k2+2*k3+k4);
end
analitik = exp(1/3*x.^3);
disp (' x y analitik')
Hasil=[x;y;analitik]'
» rk
x y analitik
Hasil =
0 1.0000 1.0000
0.1000 1.0003 1.0003
0.2000 1.0027 1.0027
0.3000 1.0090 1.0090
0.4000 1.0216 1.0216
0.5000 1.0425 1.0425
0.6000 1.0747 1.0747
0.7000 1.1211 1.1211
0.8000 1.1861 1.1861
0.9000 1.2751 1.2751
1.0000 1.3956 1.3956
6
function F = f803(x,y)
F = x*x*y;
» [x y] = ode45('f803',[0:0.1:1],1);
» [x y]
ans =
0 1.0000
0.1000 1.0003
0.2000 1.0027
0.3000 1.0090
0.4000 1.0216
0.5000 1.0425
0.6000 1.0747
0.7000 1.1211
0.8000 1.1861
0.9000 1.2751
1.0 1.3956
7
Contoh : Riggs halaman 148, problem 4.2
Sebuah reaktor semi-batch digunakan untuk reaksi fasa cair :
A ----> produk. Reaksi berlangsung dengan kecepatan r = 0,1
CA2. mula-mula reaktor hanya berisi cairan inert dengan volume 50
L. Mulai t = 0 detik, suatu larutan A dengan konsentrasi konstan 1
mol/L dan laju alir 10 L/detik dialirkan ke dalam reaktor secara
kontinu.
Pertanyaan : tentukan jumlah mol A pada saat t = 2 detik setelah
larutan A dialirkan ke dalam reaktor
Jawab :
Qo, CAo
VR(t)
VR(t)
Neraca massa untuk A :
in – consumed + generation = out + accumulation
consumed : [Link]
generation : [Link] = 0 (tidak ada generasi untuk zat A)
dn A
QoCAo - kCA2VR + 0 = 0 +
dt
CA = nA/VR
8
dn A n A2
= QoCAo - k .......................................................................... (1)
dt VR
CA = konsentrasi A dalam reaktor (mol/L)
nA = jumlah mol A dalam reaktor pada saat t (mol)
t = waktu (detik)
Ketika cairan ditambahkan ke dalam reaktor, maka volume reaktor (V R)
akan naik. ===> VR = fungsi waktu.
d ( .VR )
= Qo . ; asumsi : ρ = konstan
dt
dVR
maka : = Q ===> VR = Qo.t + C
dt
syarat batas : t = 0 -----------> VR = V0
VR = Qo.t + Vo ........................ (2)
(2) -----> (1)
dn A k .n A2
= Qo .C Ao − ; t = 0 --------> nA = 0
dt Qo .t + Vo
Dik : CAo = 1,0 mol/l
k = 0,1 l/mol.s
Qo = 10 l/s
Vo = 50 l
Maka :
dn A 0,1. n A2
= 10.1 −
dt 10.t + 50
dn A 0,1. n A2
= 10 −
dt 50 + 10 t
Metoda Euler eksplisit :
9
0,1. n A2 ( t )
nA (t + Δt) = nA(t) + t. 10 −
50 + [Link]
0,1. n A2 ,i −1
atau : nA,i = nA,i-1 + t. 10 −
50 + 10 . t i −1
t0 = 0 -------> nAo = 0 ; t = 1
0,1. n A2 , 0
t = 1 -------> nA,1 = nA,0 + t. 10 −
50 + 10 .t 0
0,1. 0 2
nA,1 = 0 + 1. 10 −
50 + 10. 0
nA,1 = 10
0,1. n A2 ,1
t = 2 -------> nA,2 = nA,1 + t. 10 −
50 + 10 .t1
0,1.10 2
nA,2 = 10 + 1. 10 −
50 + 10 .1
nA,2 = 19,833 mol
Analitik = 19,5982 --------------> dicapai jika ∆t ≈ 0,0001
10
Editor :
function F = f805(t,na)
F = 10 - (0.1*na^2)/(50 + 10*t);
Command window:
>> [t na] = ode23('f805',[0:1:2],0);
>> [t na]
ans =
0 0
1.0000 9.9424
2.0000 19.5982
11