Source code MATLAB untuk Menyelesaikan Sistem dengan 3 Persamaan (Metode Adam-Bashforth-Moulton)

%Metode Adam-Bashforth-Moulton

clc;

clear all;

h=0.01;%lebar sub interval

t=0:h:2000;%batas interval

N=length(t)-1;

S=zeros(N+1,1);

I=zeros(N+1,1);

P=zeros(N+1,1);

%nilai awal

S(1)=50;

I(1)=20;

P(1)=30;

for i=1:3

    % Runge Kutta order 4 untuk mencari titik 1 sampai 3

    k1t=f1(t(i),S(i),I(i),P(i));

    l1t=f2(t(i),S(i),I(i),P(i));

    m1t=f3(t(i),S(i),I(i),P(i));

   

    k2t=f1(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);

    l2t=f2(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);

    m2t=f3(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);

   

    k3t=f1(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);

    l3t=f2(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);

    m3t=f3(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);

   

    k4t=f1(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);

    l4t=f2(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);

    m4t=f3(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);

   

    S(i+1)=S(i)+(k1t+2*k2t+2*k3t+k4t)*(h/6);

    I(i+1)=I(i)+(l1t+2*l2t+2*l3t+l4t)*(h/6);

    P(i+1)=P(i)+(m1t+2*m2t+2*m3t+m4t)*(h/6);

end

for i=4:N

%   Adam-Bashforth-Moulton

%   Prediktor

    S(i+1)=S(i)+(h/24)*(55*f1(t(i),S(i),I(i),P(i))-59*f1(t(i-1),S(i-1),I(i-1),P(i-1))...

        +37*f1(t(i-2),S(i-2),I(i-2),P(i-2))-9*f1(t(i-3),S(i-3),I(i-3),P(i-3)));

    I(i+1)=I(i)+(h/24)*(55*f2(t(i),S(i),I(i),P(i))-59*f2(t(i-1),S(i-1),I(i-1),P(i-1))...

        +37*f2(t(i-2),S(i-2),I(i-2),P(i-2))-9*f2(t(i-3),S(i-3),I(i-3),P(i-3)));

    P(i+1)=P(i)+(h/24)*(55*f3(t(i),S(i),I(i),P(i))-59*f3(t(i-1),S(i-1),I(i-1),P(i-1))...

        +37*f3(t(i-2),S(i-2),I(i-2),P(i-2))-9*f3(t(i-3),S(i-3),I(i-3),P(i-3)));

%   Korektor

    S(i+1)=S(i)+(h/24)*(9*f1(t(i+1),S(i+1),I(i+1),P(i+1))+19*f1(t(i),S(i),I(i),P(i))...

        -5*f1(t(i-1),S(i-1),I(i-1),P(i-1))+f1(t(i-2),S(i-2),I(i-2),P(i-2)));

    I(i+1)=I(i)+(h/24)*(9*f2(t(i+1),S(i+1),I(i+1),P(i+1))+19*f2(t(i),S(i),I(i),P(i))...

        -5*f2(t(i-1),S(i-1),I(i-1),P(i-1))+f2(t(i-2),S(i-2),I(i-2),P(i-2)));

    P(i+1)=P(i)+(h/24)*(9*f3(t(i+1),S(i+1),I(i+1),P(i+1))+19*f3(t(i),S(i),I(i),P(i))...

        -5*f3(t(i-1),S(i-1),I(i-1),P(i-1))+f3(t(i-2),S(i-2),I(i-2),P(i-2)));

end

tn=t(1:N+1);

Sn=S(1:N+1)';

In=I(1:N+1)';

Pn=P(1:N+1)';

%Menampilkan hasil iterasi

fprintf('%6.3f %12.6f %18.6f %23.6f\n',[tn; Sn; In; Pn]);

%Menampilkan grafik hasil iterasi

figure(3)

plot(t,S,'r',t,I,'b',t,P,'k');

legend('S','I','P')

title('Metode Adam Bashforth Moulton Orde 4');

grid on