Senin, 19 Oktober 2015

Source Code MATLAB Metode Adam Bashforth Moulton

clc;
clear all;
%fungsi yang akan dicari penyelesaiannya
fu1=inline('9*u1+24*u2+5*cos(t)-(1/3)*sin(t)','t','u1','u2');
fu2=inline('-24*u1-51*u2-9*cos(t)+(1/3)*sin(t)','t','u1','u2');
%fungsi penyelesaian secara eksak
fu1_eks=inline('2*exp(-3*t)-exp(-39*t)+(1/3)*cos(t)','t');
fu2_eks=inline('-exp(-3*t)+2*exp(-39*t)-(1/3)*cos(t)','t');
h=0.01;%lebar sub interval
t=0:h:10;%batas interval
N=length(t)-1;
u1=zeros(N+1,1);
u2=zeros(N+1,1);
u1_eks=zeros(N+1,1);
u2_eks=zeros(N+1,1);
%nilai awal
u1(1)=4/3;
u2(1)=2/3;
u1_eks(1)=4/3;
u2_eks(1)=2/3;
error1(1)=0;
error2(1)=0;
for i=1:N
    % Runge Kutta
    k1t=fu1(t(i),u1(i),u2(i));
    l1t=fu2(t(i),u1(i),u2(i));
    k2t=fu1(t(i)+h/2,u1(i)+k1t*h/2,u2(i)+l1t*h/2);
    l2t=fu2(t(i)+h/2,u1(i)+k1t*h/2,u2(i)+l1t*h/2);
    k3t=fu1(t(i)+h/2,u1(i)+k2t*h/2,u2(i)+l2t*h/2);
    l3t=fu2(t(i)+h/2,u1(i)+k2t*h/2,u2(i)+l2t*h/2);
    k4t=fu1(t(i)+h,u1(i)+k3t*h,u2(i)+l3t*h);
    l4t=fu2(t(i)+h,u1(i)+k3t*h,u2(i)+l3t*h);
 
    u1(i+1)=u1(i)+(k1t+2*k2t+2*k3t+k4t)*(h/6);
    u2(i+1)=u2(i)+(l1t+2*l2t+2*l3t+l4t)*(h/6);
end
for i=4:N
    % Adams-Bashforth -- *predictor*
    u1(i+1) = u1(i) + h/24*(55*fu1(t(i), u1(i),u2(i)) - 59*fu1(t(i-1),...
        u1(i-1),u2(i-1))+37*fu1(t(i-2),u1(i-2),u2(i-2)) - 9*fu1(t(i-3),...
        u1(i-3),u2(i-3)));
    u2(i+1) = u2(i) + h/24*(55*fu2(t(i), u1(i),u2(i)) - 59*fu2(t(i-1),...
        u1(i-1),u2(i-1))+37*fu2(t(i-2),u1(i-2),u2(i-2)) - 9*fu2(t(i-3),...
        u1(i-3),u2(i-3)));
t(i+1) = t(i) + h;
% Adams-Moulton -- *corrector*
u1(i+1) = u1(i) + h/24*(9*fu1(t(i+1),u1(i+1),u2(i+1)) + 19*fu1(t(i),...
        u1(i),u2(i))-5*fu1(t(i-1),u1(i-1),u2(i-1)) + fu1(t(i-2),u1(i-2),...
        u2(i-2)));
    u2(i+1) = u2(i) + h/24*(9*fu2(t(i+1),u1(i+1),u2(i+1)) + 19*fu2(t(i),...
        u1(i),u2(i))-5*fu2(t(i-1),u1(i-1),u2(i-1)) + fu2(t(i-2),u1(i-2),...
        u2(i-2)));
    %solusi eksak
    u1_eks(i+1)=fu1_eks(t(i+1));
    u2_eks(i+1)=fu2_eks(t(i+1));
    %error
    error1(i+1)=abs(u1(i+1)-u1_eks(i+1));
    error2(i+1)=abs(u2(i+1)-u2_eks(i+1));
end
%menampilkan error iterasi terakhir untuk masing-masing persamaan
fprintf('nilai u1 terakhirnya %d\n',u1(N));
fprintf('nilai u2 terakhirnya %d\n',u2(N));
fprintf('Error hasil iterasi terakhir u1 adalah %d\n',error1(N));
fprintf('Error hasil iterasi terakhir u2 adalah %d\n',error2(N));
figure(1)
plot(t,u1_eks,'ko',t,u1,'y--');
legend('u1 eksak','u1');
title('Metode ABM untuk u1');xlabel('t');ylabel('u1');
figure(2)
plot(t,u2_eks,'ro',t,u2,'k--');
legend('u2 eksak','u2');
title('Metode ABM untuk u2');xlabel('t');ylabel('u2');
figure(3)
plot(t,error1,'go',t,error2,'k*');
xlabel('t');ylabel('galat u1');
legend('Galat u1','Galat u2');

Senin, 12 Oktober 2015

Source code MATLAB metode Euler Eksplisit

%Metode Euler Eksplisit
%dy/dt=f(t,y)
%f(t,y)=-(y+1)*(y+3);
clear all;
dt=0.01;
t=0:dt:2;
N=length(t)-1;
y=zeros(N+1,1);

y(1)=-2;
y_eks(1)=-2;
err(1)=0;
for n=1:N
    y(n+1)=y(n)+dt*(-(y(n)+1)*(y(n)+3));%solusi Euler eksplisit
    y_eks(n+1)=-3+2/(1+exp(-2*t(n+1)));%solusi eksak
    err(n+1)=abs(y(n+1)-y_eks(n+1));%kesalahan mutlak
end;
figure(1)
plot(t,y,'ko',t,y_eks,'r');
figure(2);%hold on
plot(t,err);

Source Code MATLAB metode Runge-Kutta Order 2 dan 4

clc;
clear all;
%fungsi yang akan dicari penyelesaiannya
fu1=inline('9*u1+24*u2+5*cos(t)-(1/3)*sin(t)','t','u1','u2');
fu2=inline('-24*u1-51*u2-9*cos(t)+(1/3)*sin(t)','t','u1','u2');
%penyelesaian secara eksak
fu1_eks=inline('2*exp(-3*t)-exp(-39*t)+(1/3)*cos(t)','t');
fu2_eks=inline('-exp(-3*t)+2*exp(-39*t)-(1/3)*cos(t)','t');
h=0.01;%lebar sub interval
t=0:h:20;%batas interval
N=length(t)-1;
u1=zeros(N+1,1);
u2=zeros(N+1,1);
u3=zeros(N+1,1);
u4=zeros(N+1,1);
u1_eks=zeros(N+1,1);
u2_eks=zeros(N+1,1);
%nilai awal
u1(1)=4/3;
u2(1)=2/3;
u3(1)=4/3;
u4(1)=2/3;
u1_eks(1)=4/3;
u2_eks(1)=2/3;
u3_eks(1)=4/3;
u4_eks(1)=2/3;
error1(1)=0;
error2(1)=0;
for i=1:N
    %Metode Runge Kutta Order 2
    k1=h*fu1(t(i),u1(i),u2(i));
    l1=h*fu2(t(i),u1(i),u2(i));
    k2=h*fu1(t(i)+h,u1(i)+k1,u2(i)+l1);
    l2=h*fu2(t(i)+h,u1(i)+k1,u2(i)+l1);
    %solusi numerik
    u1(i+1)=u1(i)+(k1+k2)/2;
    u2(i+1)=u2(i)+(l1+l2)/2;
    %solusi eksak
    u1_eks(i+1)=fu1_eks(t(i+1));
    u2_eks(i+1)=fu2_eks(t(i+1));
    %error
    error1(i+1)=abs(u1(i+1)-u1_eks(i+1));
    error2(i+1)=abs(u2(i+1)-u2_eks(i+1));
    %Metode Runge Kutta Order 4
    k1t=fu1(t(i),u3(i),u4(i));
    l1t=fu2(t(i),u3(i),u4(i));
    k2t=fu1(t(i)+h/2,u3(i)+k1t*h/2,u4(i)+l1t*h/2);
    l2t=fu2(t(i)+h/2,u3(i)+k1t*h/2,u4(i)+l1t*h/2);
    k3t=fu1(t(i)+h/2,u3(i)+k2t*h/2,u4(i)+l2t*h/2);
    l3t=fu2(t(i)+h/2,u3(i)+k2t*h/2,u4(i)+l2t*h/2);
    k4t=fu1(t(i)+h,u3(i)+k3t*h,u4(i)+l3t*h);
    l4t=fu2(t(i)+h,u3(i)+k3t*h,u4(i)+l3t*h);
    %solusi numerik
    u3(i+1)=u3(i)+(k1t+2*k2t+2*k3t+k4t)*(h/6);
    u4(i+1)=u4(i)+(l1t+2*l2t+2*l3t+l4t)*(h/6);
    %solusi eksak
    u3_eks(i+1)=fu1_eks(t(i+1));
    u4_eks(i+1)=fu2_eks(t(i+1));
    %error
    error3(i+1)=abs(u3(i+1)-u3_eks(i+1));
    error4(i+1)=abs(u4(i+1)-u4_eks(i+1));
end
%menampilkan error iterasi terakhir untuk masing-masing persamaan
fprintf('Error hasil iterasi terakhir u1 metode RK2 adalah %d\n',error1(N));
fprintf('Error hasil iterasi terakhir u1 metode RK4 adalah %d\n',error3(N));
fprintf('Error hasil iterasi terakhir u2 metode RK2 adalah %d\n',error2(N));
fprintf('Error hasil iterasi terakhir u2 metode RK2 adalah %d\n',error4(N));
figure(1)
plot(t,u1_eks,'ko',t,u1,'b*',t,u3,'g--');
legend('u1 eksak','u1 metode RK2','u1 metode RK4');
title('Metode Runge-Kutta untuk u1');xlabel('t');ylabel('u1');
figure(2)
plot(t,u2_eks,'ro',t,u2,'k*',t,u4,'y--');
legend('u2 eksak','u2 metode RK2','u2 metode RK4');
title('Metode Runga-Kutta untuk u2');xlabel('t');ylabel('u2');
figure(3)
plot(t,error1,'go',t,error3,'k*');
xlabel('t');ylabel('galat u1');
legend('Galat u1 metode RK2','Galat u1 metode RK4');
figure(4)
plot(t,error2,'ro',t,error4,'b*');
xlabel('t');ylabel(' galat u2');
legend('Galat u2 metode RK2','Galat u2 metode RK4');