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');