NAMA : AGUS ARDIANSYAH
NIM : 16524043
1. Metode euler
a. Percobaan 1
function [T,Y] = myeuler(f,tspan,y0,n)
a = tspan(1); b = tspan(2);
h = (b-a)/n;
t = a; T=a;
y = y0;
Y = y0';
for i = 1: n
y = y+h*f(t,y);
t = a+i*h;
T = [T;t];
Y = [Y;y'];
end
end
g = inline('[y(2);-y(1)]','t','y');
N1 = 5;
[T Y] = myeuler(g,[0 4*pi],[1;0],N1);
y1 = Y(:,1);
y2 = Y(:,2);
figure
subplot 121
plot(y1,y2);
N1 = 20
N1=100
N1 = 200
N1 = 1000
N1 = 10000
Analisa dan kesimpulan : Pada program ini yang sangat mempengaruhi hasil akhir adalah
nilai inputan dari N1. Sehingga untuk menghasilkan gambar circle yang lebih baik maka
nilai N1 harus tinggi. Bisa dilihat dari beberapa percobaan diatas, semakin tinggi nilai
N1 maka gambar yang dihasilkan lebih baik dan mendekati circle. Jadi, percobaan ini
adalah termasuk unit circle.
b. Percobaan 2
function [T,Y] = mymodeuler (f,tspan,y0,n)
a = tspan (1); b = tspan (2);
h = (b - a )/ n;
t = a; T = a ;
y = y0 ;
Y = y0';
for i = 1: n
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
y = y +0.5*(k1+k2);
t = a+i*h;
T = [T;t];
Y = [Y;y'];
end
end
g = inline('[y(2);-y(1)]','t','y');
N1 = 25;
[T Y] = myeuler(g,[0 4*pi],[1;0],N1);
y1 = Y(:,1);
y2 = Y(:,2);
figure
subplot 121
plot(y1,y2);
N2 = 25;
[Tm Ym] = mymodeuler(g,[0 4*pi],[1;0],N2);
ym1 = Ym(:,1);
ym2 = Ym(:,2);
subplot 122
plot(ym1,ym2)
N2=20
N2=100
N2= 1000
N2=10000
Analisa dan kesimpulan : dengan menggunakan percobaan 2 maka hal yang paling
diuntungkan adalah untuk mendapatkan hasil yang diinginkan tidak perlu melakukan
banyak pendeketan. Jadi, bisa dikatakan percobaan 2 ini memiliki error yang lebih sedikit
deibandingkan dengan percobaan 1.
2. Metode Runge-kutta
a. Fungsi euler
function [T,Y] = mymodeuler (f,tspan,y0,n)
a = tspan (1); b = tspan (2);
h = (b - a )/ n;
t = a; T = a ;
y = y0 ;
Y = y0';
for i = 1: n
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
y = y +0.5*(k1+k2);
t = a+i*h;
T = [T;t];
Y = [Y;y'];
end
end
clear all; close all; clc;
f = inline('[y(2);-.1*y(2)-sin(y(1))+sin(t)]','t','y');
N = 20;
[T Y] = myeuler(f,[0 20],[1;-1.5],N);
size(T);
size(Y);
theta = Y(:,1);
figure;
plot(T,theta); hold on
[Te Ye]= ode45(f,[0 20],[1;-1.5]);
theta = Ye(:,1);
plot(Te,theta,'r')
legend('Euler','Runge-Kutta')
N= 100
N=200
N=1000
N=10000
b. Fungsi mymodeuler
function [T,Y] = mymodeuler (f,tspan,y0,n)
a = tspan (1); b = tspan (2);
h = (b - a )/ n;
t = a; T = a ;
y = y0 ;
Y = y0';
for i = 1: n
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
y = y +0.5*(k1+k2);
t = a+i*h;
T = [T;t];
Y = [Y;y'];
end
end
clear all; close all; clc;
f = inline('[y(2);-.1*y(2)-sin(y(1))+sin(t)]','t','y');
N = 20;
[T Y] = mymodeuler(f,[0 20],[1;-1.5],N);
size(T);
size(Y);
theta = Y(:,1);
figure;
plot(T,theta); hold on
[Te Ye]= ode45(f,[0 20],[1;-1.5]);
theta = Ye(:,1);
plot(Te,theta,'r')
legend('Euler','Runge-Kutta')
N=100
N=200
N=1000
Analisa dan kesimpulan : pada percobaan ini terdapat 2 fungsi yaitu fungsi euler dan
fungsi mymodeuler, saat melakukan percobaan dan mendapatkan hasil maka yang
mempengaruhi adalah inputan nilai N nya. Sama hal nya seperti sebelumnya, semaknin
nilai N dinakaikan dari sebelumnya maka eror yang dihasilkan akan lebih sedikit. Jika
dibandingkan keduanya maka fungsi yang lebih baik adalah mymodeuler karena dengan
menggunkan mymodeuler hasil yang didapat lebih akurat. Saya ambil contoh pada fungsi
euler grafik yang dihasilkan sehingga grafik sama dengan rungge kutta adalah pada nilai
N=1000. Jika pada mymodeuler hanya butuh masukan nilai N=200 untuk mendapatkan
hasil yang sama. Jadi, kesimpulannya adalah saat kita bandingkan kedua fungsi maka
fungsi yang lebih baik dalam kinerja adalah fungsi mymodeuler.
3. Bentuk lain penulisan fungsi (contoh pada kasus Lorentz attractor)
clear all; close all; clc;
f=@(t,x,beta,rho,sigma)[-beta*x(1)+x(2)*x(3);...
-rho*x(2)+rho*x(3); -x(1)*x(2)+sigma*x(2)-x(3)];
t_final=100; x0=[0;0;1e-10];
b1=8/3; r1=10; s1=28;
[t,x]=ode45(f,[0,t_final],x0,[],b1,r1,s1);
figure
plot(t,x);
figure;
plot3(x(:,1),x(:,2),x(:,3));
axis([10 42 -20 20 -20 28]);
comet3(x(:,1),x(:,2),x(:,3));
Analisa dan kesimpulan : pada percobaan ini menghasilkan gambar 3D seperti halnya ada
2 linggkaran yang berdekatan tetapi lingkaran yang satu tampak berbaring dan satunya
tampak berdiri. Hal ini dipengaruhi oleh nilai inputan dari beta, rho, dan sigma.
Pengaruhnya adalah saat kita mengganti nilai beta dan rho lebih tepatnya mengecilkan
nilainya dari sebelumnya maka akan menghasilkan gambar yang berbeda dan
lingkarannya menjadi 1 saja. Contoh: b1=1; r1=2; s1=10
4. Tugas soal
a. Tugas 1
clear all; close all; clc;
f=@(t,x,beta,rho,sigma)[-beta*x(1)+x(2)*x(3);...
-rho*x(2)+rho*x(3); -x(1)*x(2)+sigma*x(2)-x(3)];
t_final=100; x0=[0;0;1e-10];
b1=2; r1=5; s1=20;
[t,x]=ode45(f,[0,t_final],x0,[],b1,r1,s1);
figure
plot(t,x);
figure;
plot3(x(:,1),x(:,2),x(:,3));
axis([10 42 -20 20 -20 28]);
comet3(x(:,1),x(:,2),x(:,3));
b. Tugas 2
M=1
clear all; close all; clc;
f=@(t,x,m)[x(2);...
-m*(x(1)^2-1)*x(2)-x(1)];
t_final=20; x0=[-0.2 -0.7];
m=1;
[t,x]=ode45(f,[0,t_final],x0,[],m);
figure
plot(t,x);
figure;
M=2
clear all; close all; clc;
f=@(t,x,m)[x(2);...
-m*(x(1)^2-1)*x(2)-x(1)];
t_final=20; x0=[-0.2 -0.7];
m=2;
[t,x]=ode45(f,[0,t_final],x0,[],m);
figure
plot(t,x);
figure;
Analisa dan kesimpulan : dari apa yang telah dilakukan sebelumnya terdapat 2 input nilai
m yang berbeda, yaitu m=1 dan m=2. Sebelumnya, m adalah mewakili nilai mikro. Dari
hasil yang didapat dari nilai input 1 dan 2 terdapat perbedaan yaitu error yang dihasilkan
lebih baik pada m=2. Jika dilihat dari kasat mata perbedaannya hanya sedikit saja.
Walaupun demikian pada hasil gambar di kedua inputan tersebut masih memiliki error
yang cukup besar. Jadi, untuk membuktikan bahwasanya nilai m yang mempengaruhi hasil
akhirnya maka saya melakukan tambahan percobaan dengan menggantikan nilai m, lebih
tepatnya menaikkan nilainya menjadi 1000.
Bisa dilihat gambar menunjukan error yang jauh lebih kecil disbanding sebelumnya. Oleh
sebab itu, yang mempengaruhi hasil akhir adalah nilai mikronya. Jadi, kesimpulannya
untuk mendapatkan hasil yang errornya kecil maka nilai mikro harus dinaikkan dari
sebelumnya.
