Please correct my MATLAB code about Seventh Order Adam-Bashforth-Moulton for solving ODES (Edit: With ANOTHER EXAMPLE)

조회 수: 1 (최근 30일)
Fahmy Shandy
Fahmy Shandy 2019년 12월 13일
댓글: Looky 2019년 12월 13일
clear all;
close all;
clc;
disp('=======================================');
disp(' ADAM-BASHFORTH-MOULTON SEVEN STEPS ');
disp(' WITH SEVENTH ORDER RUNGE-KUTTA ');
disp(' by Shandy Verdyo ');
disp('=======================================');
% ORDINARY DIFFERENTIAL EQUATIONS SYSTEM
fS = @(t,S,E,I,R)(-(0.5*I + 0.25)*S);
fE = @(t,S,E,I,R)(0.5*S*I - E);
fI = @(t,S,E,I,R)(0.75*E - I);
fR = @(t,S,E,I,R)(0.75*I - 0.25*R);
t0 = input('Input t (Input 0.25)= ');
S0 = input('Input S0 (Input 1)= ');
E0 = input('Input E0 (Input 0)= ');
I0 = input('Input I0 (Input 1)= ');
R0 = input('Input R0 (Input 0)= ');
h = 0.01;
n = (t0)/h;
t(1) = t0;
S(1) = S0;
E(1) = E0;
I(1) = I0;
R(1) = R0;
disp('_________________________________________________________________');
fprintf(' t S E I R \n');
disp('_________________________________________________________________');
fprintf(' %5.4f %11.8f %11.8f %11.8f %11.8f\n',t(1),S(1),...
E(1),I(1),R(1));
for i = 1:6 %Seveth order Runge-Kutta for finding seven initial value
t(i+1) = t(i)+h;
k1S = fS(t(i),S(i),E(i),I(i),R(i));
k1E = fE(t(i),S(i),E(i),I(i),R(i));
k1I = fI(t(i),S(i),E(i),I(i),R(i));
k1R = fR(t(i),S(i),E(i),I(i),R(i));
k2S = fS(t(i)+h/12, S(i)+h*k1S/12, E(i)+h*k1E/12, I(i)+h*k1I/12,...
R(i)+h*k1R/12);
k2E = fE(t(i)+h/12, S(i)+h*k1S/12, E(i)+h*k1E/12, I(i)+h*k1I/12,...
R(i)+h*k1R/12);
k2I = fI(t(i)+h/12, S(i)+h*k1S/12, E(i)+h*k1E/12, I(i)+h*k1I/12,...
R(i)+h*k1R/12);
k2R = fR(t(i)+h/12, S(i)+h*k1S/12, E(i)+h*k1E/12, I(i)+h*k1I/12,...
R(i)+h*k1R/12);
k3S = fS(t(i)+h/6, S(i)+h/6*(k2S), E(i)+h/6*(k2E), I(i)+h/6*(k2I),...
R(i)+h/6*(k2R));
k3E = fE(t(i)+h/6, S(i)+h/6*(k2S), E(i)+h/6*(k2E), I(i)+h/6*(k2I),...
R(i)+h/6*(k2R));
k3I = fI(t(i)+h/6, S(i)+h/6*(k2S), E(i)+h/6*(k2E), I(i)+h/6*(k2I),...
R(i)+h/6*(k2R));
k3R = fR(t(i)+h/6, S(i)+h/6*(k2S), E(i)+h/6*(k2E), I(i)+h/6*(k2I),...
R(i)+h/6*(k2R));
k4S = fS(t(i)+h/4, S(i)+h/16*(k1S+3*k3S), E(i)+h/16*(k1E+3*k3E),...
I(i)+h/16*(k1I+3*k3I), R(i)+h/16*(k1R+3*k3R));
k4E = fE(t(i)+h/4, S(i)+h/16*(k1S+3*k3S), E(i)+h/16*(k1E+3*k3E),...
I(i)+h/16*(k1I+3*k3I), R(i)+h/16*(k1R+3*k3R));
k4I = fI(t(i)+h/4, S(i)+h/16*(k1S+3*k3S), E(i)+h/16*(k1E+3*k3E),...
I(i)+h/16*(k1I+3*k3I), R(i)+h/16*(k1R+3*k3R));
k4R = fR(t(i)+h/4, S(i)+h/16*(k1S+3*k3S), E(i)+h/16*(k1E+3*k3E),...
I(i)+h/16*(k1I+3*k3I), R(i)+h/16*(k1R+3*k3R));
k5S = fS(t(i)+3*h/4, S(i)+h/16*(21*k1S - 81*k3S + 72*k4S), E(i)+...
h/16*(21*k1E - 81*k3E + 72*k4E),I(i)+h/16*(21*k1I - 81*k3I +...
72*k4I), R(i)+h/16*(21*k1R - 81*k3R + 72*k4R));
k5E = fE(t(i)+3*h/4, S(i)+h/16*(21*k1S - 81*k3S + 72*k4S), E(i)+...
h/16*(21*k1E - 81*k3E + 72*k4E),I(i)+h/16*(21*k1I - 81*k3I +...
72*k4I), R(i)+h/16*(21*k1R - 81*k3R + 72*k4R));
k5I = fI(t(i)+3*h/4, S(i)+h/16*(21*k1S - 81*k3S + 72*k4S), E(i)+...
h/16*(21*k1E - 81*k3E + 72*k4E),I(i)+h/16*(21*k1I - 81*k3I +...
72*k4I), R(i)+h/16*(21*k1R - 81*k3R + 72*k4R));
k5R = fR(t(i)+3*h/4, S(i)+h/16*(21*k1S - 81*k3S + 72*k4S), E(i)+...
h/16*(21*k1E - 81*k3E + 72*k4E),I(i)+h/16*(21*k1I - 81*k3I +...
72*k4I), R(i)+h/16*(21*k1R - 81*k3R + 72*k4R));
k6S = fS(t(i)+16*h/17, S(i)+h/250563*(1344688*k1S...
- 5127552*k3S + 4096896*k4S - 78208*k5S), E(i)+h/250563*...
(1344688*k1E - 5127552*k3E + 4096896*k4E - 78208*k5E), I(i)+...
h/250563*(1344688*k1I - 5127552*k3I + 4096896*k4I - 78208*k5I),...
R(i)+h/250563*(1344688*k1R - 5127552*k3R + 4096896*k4R -...
78208*k5R));
k6E = fE(t(i)+16*h/17, S(i)+h/250563*(1344688*k1S...
- 5127552*k3S + 4096896*k4S - 78208*k5S), E(i)+h/250563*...
(1344688*k1E - 5127552*k3E + 4096896*k4E - 78208*k5E), I(i)+...
h/250563*(1344688*k1I - 5127552*k3I + 4096896*k4I - 78208*k5I),...
R(i)+h/250563*(1344688*k1R - 5127552*k3R + 4096896*k4R -...
78208*k5R));
k6I = fI(t(i)+16*h/17, S(i)+h/250563*(1344688*k1S...
- 5127552*k3S + 4096896*k4S - 78208*k5S), E(i)+h/250563*...
(1344688*k1E - 5127552*k3E + 4096896*k4E - 78208*k5E), I(i)+...
h/250563*(1344688*k1I - 5127552*k3I + 4096896*k4I - 78208*k5I),...
R(i)+h/250563*(1344688*k1R - 5127552*k3R + 4096896*k4R -...
78208*k5R));
k6R = fR(t(i)+16*h/17, S(i)+h/250563*(1344688*k1S...
- 5127552*k3S + 4096896*k4S - 78208*k5S), E(i)+h/250563*...
(1344688*k1E - 5127552*k3E + 4096896*k4E - 78208*k5E), I(i)+...
h/250563*(1344688*k1I - 5127552*k3I + 4096896*k4I - 78208*k5I),...
R(i)+h/250563*(1344688*k1R - 5127552*k3R + 4096896*k4R -...
78208*k5R));
k7S = fS(t(i)+h/2, S(i)+h/234624*( -341549*k1S + 1407744*k3S...
- 1018368*k4S + 84224*k5S - 14739*k6S), E(i)+h/234624*...
( -341549*k1E + 1407744*k3E - 1018368*k4E + 84224*k5E -...
14739*k6E), I(i)+h/234624*( -341549*k1I + 1407744*k3I -...
1018368*k4I + 84224*k5I - 14739*k6I), R(i)+h/234624*...
( -341549*k1R + 1407744*k3R - 1018368*k4R + 84224*k5R -...
14739*k6R));
k7E = fE(t(i)+h/2, S(i)+h/234624*( -341549*k1S + 1407744*k3S...
- 1018368*k4S + 84224*k5S - 14739*k6S), E(i)+h/234624*...
( -341549*k1E + 1407744*k3E - 1018368*k4E + 84224*k5E -...
14739*k6E), I(i)+h/234624*( -341549*k1I + 1407744*k3I -...
1018368*k4I + 84224*k5I - 14739*k6I), R(i)+h/234624*...
( -341549*k1R + 1407744*k3R - 1018368*k4R + 84224*k5R -...
14739*k6R));
k7I = fI(t(i)+h/2, S(i)+h/234624*( -341549*k1S + 1407744*k3S...
- 1018368*k4S + 84224*k5S - 14739*k6S), E(i)+h/234624*...
( -341549*k1E + 1407744*k3E - 1018368*k4E + 84224*k5E -...
14739*k6E), I(i)+h/234624*( -341549*k1I + 1407744*k3I -...
1018368*k4I + 84224*k5I - 14739*k6I), R(i)+h/234624*...
( -341549*k1R + 1407744*k3R - 1018368*k4R + 84224*k5R -...
14739*k6R));
k7R = fR(t(i)+h/2, S(i)+h/234624*( -341549*k1S + 1407744*k3S...
- 1018368*k4S + 84224*k5S - 14739*k6S), E(i)+h/234624*...
( -341549*k1E + 1407744*k3E - 1018368*k4E + 84224*k5E -...
14739*k6E), I(i)+h/234624*( -341549*k1I + 1407744*k3I -...
1018368*k4I + 84224*k5I - 14739*k6I), R(i)+h/234624*...
( -341549*k1R + 1407744*k3R - 1018368*k4R + 84224*k5R -...
14739*k6R));
k8S = fS(t(i)+h, S(i)+h/136864*( -381875*k1S + 1642368*k3S...
- 1327872*k4S + 72192*k5S + 14739*k6S + 117312*k7S), E(i)+...
h/136864*( -381875*k1E + 1642368*k3E - 1327872*k4E + 72192*k5E +...
14739*k6E + 117312*k7E), I(i)+h/136864*( -381875*k1I +...
1642368*k3I - 1327872*k4I + 72192*k5I + 14739*k6I + 117312*k7I),...
R(i)+h/136864*( -381875*k1R + 1642368*k3R - 1327872*k4R +...
72192*k5R + 14739*k6R + 117312*k7R));
k8E = fE(t(i)+h, S(i)+h/136864*( -381875*k1S + 1642368*k3S...
- 1327872*k4S + 72192*k5S + 14739*k6S + 117312*k7S), E(i)+...
h/136864*( -381875*k1E + 1642368*k3E - 1327872*k4E + 72192*k5E +...
14739*k6E + 117312*k7E), I(i)+h/136864*( -381875*k1I +...
1642368*k3I - 1327872*k4I + 72192*k5I + 14739*k6I + 117312*k7I),...
R(i)+h/136864*( -381875*k1R + 1642368*k3R - 1327872*k4R +...
72192*k5R + 14739*k6R + 117312*k7R));
k8I = fI(t(i)+h, S(i)+h/136864*( -381875*k1S + 1642368*k3S...
- 1327872*k4S + 72192*k5S + 14739*k6S + 117312*k7S), E(i)+...
h/136864*( -381875*k1E + 1642368*k3E - 1327872*k4E + 72192*k5E +...
14739*k6E + 117312*k7E), I(i)+h/136864*( -381875*k1I +...
1642368*k3I - 1327872*k4I + 72192*k5I + 14739*k6I + 117312*k7I),...
R(i)+h/136864*( -381875*k1R + 1642368*k3R - 1327872*k4R +...
72192*k5R + 14739*k6R + 117312*k7R));
k8R = fR(t(i)+h, S(i)+h/136864*( -381875*k1S + 1642368*k3S...
- 1327872*k4S + 72192*k5S + 14739*k6S + 117312*k7S), E(i)+...
h/136864*( -381875*k1E + 1642368*k3E - 1327872*k4E + 72192*k5E +...
14739*k6E + 117312*k7E), I(i)+h/136864*( -381875*k1I +...
1642368*k3I - 1327872*k4I + 72192*k5I + 14739*k6I + 117312*k7I),...
R(i)+h/136864*( -381875*k1R + 1642368*k3R - 1327872*k4R +...
72192*k5R + 14739*k6R + 117312*k7R));
k9S = fS(t(i)+2*h/3, S(i)+h/16755336*( -2070757*k1S + 9929088*k3S...
+ 584064*k4S + 3023488*k5S - 447083*k6S + 151424*k7S), E(i)+...
h/16755336*( -2070757*k1E + 9929088*k3E + 584064*k4E +...
3023488*k5E - 447083*k6E + 151424*k7E), I(i)+h/16755336*...
( -2070757*k1I + 9929088*k3I + 584064*k4I + 3023488*k5I -...
447083*k6I + 151424*k7I), R(i)+h/16755336*( -2070757*k1R +...
9929088*k3R + 584064*k4R + 3023488*k5R - 447083*k6R + 151424*k7R));
k9E = fE(t(i)+2*h/3, S(i)+h/16755336*( -2070757*k1S + 9929088*k3S...
+ 584064*k4S + 3023488*k5S - 447083*k6S + 151424*k7S), E(i)+...
h/16755336*( -2070757*k1E + 9929088*k3E + 584064*k4E +...
3023488*k5E - 447083*k6E + 151424*k7E), I(i)+h/16755336*...
( -2070757*k1I + 9929088*k3I + 584064*k4I + 3023488*k5I -...
447083*k6I + 151424*k7I), R(i)+h/16755336*( -2070757*k1R +...
9929088*k3R + 584064*k4R + 3023488*k5R - 447083*k6R + 151424*k7R));
k9I = fI(t(i)+2*h/3, S(i)+h/16755336*( -2070757*k1S + 9929088*k3S...
+ 584064*k4S + 3023488*k5S - 447083*k6S + 151424*k7S), E(i)+...
h/16755336*( -2070757*k1E + 9929088*k3E + 584064*k4E +...
3023488*k5E - 447083*k6E + 151424*k7E), I(i)+h/16755336*...
( -2070757*k1I + 9929088*k3I + 584064*k4I + 3023488*k5I -...
447083*k6I + 151424*k7I), R(i)+h/16755336*( -2070757*k1R +...
9929088*k3R + 584064*k4R + 3023488*k5R - 447083*k6R + 151424*k7R));
k9R = fR(t(i)+2*h/3, S(i)+h/16755336*( -2070757*k1S + 9929088*k3S...
+ 584064*k4S + 3023488*k5S - 447083*k6S + 151424*k7S), E(i)+...
h/16755336*( -2070757*k1E + 9929088*k3E + 584064*k4E +...
3023488*k5E - 447083*k6E + 151424*k7E), I(i)+h/16755336*...
( -2070757*k1I + 9929088*k3I + 584064*k4I + 3023488*k5I -...
447083*k6I + 151424*k7I), R(i)+h/16755336*( -2070757*k1R +...
9929088*k3R + 584064*k4R + 3023488*k5R - 447083*k6R + 151424*k7R));
k10S= fS(t(i)+h, S(i)+h/10743824*( 130521209*k1S - 499279872*k3S -...
391267968*k4S + 13012608*k5S - 3522621*k6S + 9033024*k7S...
- 30288492*k9S), E(i)+h/10743824*( 130521209*k1E - 499279872*k3E -...
391267968*k4E + 13012608*k5E - 3522621*k6E + 9033024*k7E...
- 30288492*k9E), I(i)+h/10743824*( 130521209*k1I - 499279872*k3I -...
391267968*k4I + 13012608*k5I - 3522621*k6I + 9033024*k7I...
- 30288492*k9I), R(i)+h/10743824*( 130521209*k1R - 499279872*k3R -...
391267968*k4R + 13012608*k5R - 3522621*k6R + 9033024*k7R...
- 30288492*k9R));
k10E= fS(t(i)+h, S(i)+h/10743824*( 130521209*k1S - 499279872*k3S -...
391267968*k4S + 13012608*k5S - 3522621*k6S + 9033024*k7S...
- 30288492*k9S), E(i)+h/10743824*( 130521209*k1E - 499279872*k3E -...
391267968*k4E + 13012608*k5E - 3522621*k6E + 9033024*k7E...
- 30288492*k9E), I(i)+h/10743824*( 130521209*k1I - 499279872*k3I -...
391267968*k4I + 13012608*k5I - 3522621*k6I + 9033024*k7I...
- 30288492*k9I), R(i)+h/10743824*( 130521209*k1R - 499279872*k3R -...
391267968*k4R + 13012608*k5R - 3522621*k6R + 9033024*k7R...
- 30288492*k9R));
k10I= fS(t(i)+h, S(i)+h/10743824*( 130521209*k1S - 499279872*k3S -...
391267968*k4S + 13012608*k5S - 3522621*k6S + 9033024*k7S...
- 30288492*k9S), E(i)+h/10743824*( 130521209*k1E - 499279872*k3E -...
391267968*k4E + 13012608*k5E - 3522621*k6E + 9033024*k7E...
- 30288492*k9E), I(i)+h/10743824*( 130521209*k1I - 499279872*k3I -...
391267968*k4I + 13012608*k5I - 3522621*k6I + 9033024*k7I...
- 30288492*k9I), R(i)+h/10743824*( 130521209*k1R - 499279872*k3R -...
391267968*k4R + 13012608*k5R - 3522621*k6R + 9033024*k7R...
- 30288492*k9R));
k10R= fS(t(i)+h, S(i)+h/10743824*( 130521209*k1S - 499279872*k3S -...
391267968*k4S + 13012608*k5S - 3522621*k6S + 9033024*k7S...
- 30288492*k9S), E(i)+h/10743824*( 130521209*k1E - 499279872*k3E -...
391267968*k4E + 13012608*k5E - 3522621*k6E + 9033024*k7E...
- 30288492*k9E), I(i)+h/10743824*( 130521209*k1I - 499279872*k3I -...
391267968*k4I + 13012608*k5I - 3522621*k6I + 9033024*k7I...
- 30288492*k9I), R(i)+h/10743824*( 130521209*k1R - 499279872*k3R -...
391267968*k4R + 13012608*k5R - 3522621*k6R + 9033024*k7R...
- 30288492*k9R));
S(i+1) = S(i)+h/90*(7*k1S+32*k4S+32*k5S+12*k7S+7*k8S);
E(i+1) = S(i)+h/90*(7*k1E+32*k4E+32*k5E+12*k7E+7*k8E);
I(i+1) = S(i)+h/90*(7*k1I+32*k4I+32*k5I+12*k7I+7*k8I);
R(i+1) = S(i)+h/90*(7*k1R+32*k4R+32*k5R+12*k7R+7*k8R);
fprintf('%4.f %5.4f %6.4f %6.4f %6.4f %6.4f\n',...
i,t(i),S(i),E(i),I(i),R(i));
plot(t,S,'r')
hold on
plot(t,E,'y')
hold on
plot(t,I,'g')
hold on
plot(t,R,'k')
hold on
title('Adam-Bashforth-Method'); grid on
xlabel('time(10 years)','Fontsize',10)
ylabel('the number of computers ','Fontsize',10)
legend('Susceptible','Exposed','Infectiuos','Recovered')
end
for i = 7:n %Adam-Bashforth and Adam-Moulton for finding 8th value until n
tt0 = t0+i*h;
predS = S(7) + h*(198721*fS(t(7),S(7)) - 447288*fS(t(6),S(6)) + 705549*...
fS(t(5),S(5)) - 688256*fS(t(4),S(4)) + 407139*fS(t(3),S(3)) -...
134472*fS(t(2),S(2)) + 19087*fS(t(1),S(1)))/60480;
predE = E(7) + h*(198721*fE(t(7),E(7)) - 447288*fE(t(6),E(6)) + 705549*...
fE(t(5),E(5)) - 688256*fE(t(4),E(4)) + 407139*fE(t(3),E(3)) -...
134472*fE(t(2),E(2)) + 19087*fE(t(1),E(1)))/60480;
predI = I(7) + h*(198721*fI(t(7),I(7)) - 447288*fI(t(6),I(6)) + 705549*...
fI(t(5),I(5)) - 688256*fI(t(4),I(4)) + 407139*fI(t(3),I(3)) -...
134472*fI(t(2),I(2)) + 19087*fI(t(1),I(1)))/60480;
predR = R(7) + h*(198721*fR(t(7),R(7)) - 447288*fR(t(6),R(6)) + 705549*...
fR(t(5),R(5)) - 688256*fR(t(4),R(4)) + 407139*fR(t(3),R(3)) -...
134472*fR(t(2),R(2)) + 19087*fR(t(1),R(1)))/60480;
predS = S(7) + h*(19087*fS(tt0,predS) + 65112*fS(t(7),S(7)) - 46461*...
fS(t(6),S(6)) + 37504*fS(t(5),S(5)) - 20211*fS(t(4),S(4)) + 6312*...
fS(t(3),S(3)) - 863*fS(t(2),S(2)))/60480;
predE = E(7) + h*(19087*fE(tt0,predE) + 65112*fE(t(7),E(7)) - 46461*...
fE(t(6),E(6)) + 37504*fE(t(5),E(5)) - 20211*fE(t(4),E(4)) + 6312*...
fE(t(3),E(3)) - 863*fE(t(2),E(2)))/60480;
predI = I(7) + h*(19087*fI(tt0,predI) + 65112*fI(t(7),I(7)) - 46461*...
fI(t(6),I(6)) + 37504*fI(t(5),I(5)) - 20211*fI(t(4),I(4)) + 6312*...
fI(t(3),I(3)) - 863*fI(t(2),I(2)))/60480;
predR = R(7) + h*(19087*fR(tt0,predR) + 65112*fR(t(7),R(7)) - 46461*...
fR(t(6),R(6)) + 37504*fR(t(5),R(5)) - 20211*fR(t(4),R(4)) + 6312*...
fR(t(3),R(3)) - 863*fR(t(2),R(2)))/60480;
fprintf('%4.f %5.4f %6.4f %6.4f %6.4f %6.4f\n',...
i,tt0,predS,predE,predI,predR);
for j = 1:6
t(j) = t(j+1);
S(j) = S(j+1);
E(j) = E(j+1);
I(j) = I(j+1);
R(j) = R(j+1);
end
t(7)=tt0;
S(7)=predS;
E(7)=predE;
I(7)=predI;
R(7)=predR;
plot(t,S,'r')
hold on
plot(t,E,'y')
hold on
plot(t,I,'g')
hold on
plot(t,R,'k')
title('Adam-Bashforth-Method'); grid on
xlabel('time(10 years)','Fontsize',10)
ylabel('the number of computers ','Fontsize',10)
legend('Susceptible','Exposed','Infectiuos','Recovered')
end
hold off
disp(['the number of iteration = ',num2str(i)]);
% fprintf('The result of yS with 7th order Adam-Bashforth-Moulton= %11.10f\n'...
% ,predS);
disp(['The result of yS with 7th order Adam-Bashforth-Moulton= ', num2str(predS)]);
% fprintf('The result of yE with 7th order Adam-Bashforth-Moulton= %11.10f\n'...
% ,predE);
disp(['The result of yE with 7th order Adam-Bashforth-Moulton= ', num2str(predE)]);
% fprintf('The result of yI with 7th order Adam-Bashforth-Moulton= %11.10f\n'...
% ,predI);
disp(['The result of yI with 7th order Adam-Bashforth-Moulton= ', num2str(predI)]);
% fprintf('Hasil yR dengan Adam-Bashforth-Moulton orde 7= %11.10f\n'...
% ,predR);
disp(['The result of yR with 7th order Adam-Bashforth-Moulton= ', num2str(predR)]);
I know this is very long script that might be bored for you. But a lot of time i spend for it, and it still ERROR. That's what made me cry. Please do a correction if there are several mistakes (may be). But i don't know, earlier, i'm running my code and reduce the ODEs to the ODE (just one equation) and it works. And trust me, there is no mistake about The Runge-Kutta, and ADM formula because i checked it several times and compared it with the formula that i got on the internet.
So, don't worry about RK, and ADM. Please just focused on the error.
I'm very very very grateful if you want to help me to fix this until it can shows the output.
Anyway' i've mentioned about the input above (in the script) so don't worry abou it.
For your convenient, i'll write it down my ODEs and what are its initial conditions, the equations, the input, and the ouput.
My initial conditions is :
The ODEs is :
When you input The output must be :
I hope you want to help me. I really expected an answer for my question. And i will thanks for you later. :)
clear all;
close all;
clc;
% PERSAMAAN DIFERENSIAL BIASA
f = @(x,y)(y-x^2+1);
x0 = input('Masukkan nilai x0 (Input 0)= ');
y0 = input('Masukkan nilai y0 (Input 0.5)= ');
x1 = input('Masukkan nilai x1 yang dicari (Input 2)= ');
h = input('Masukkan Ukuran Langkah <= 0.01 (Input 0.01): ');
n = (x1-x0)/h;
x(1)=x0;
y(1)=y0;
disp('____________________________');
fprintf('x y(Runge-Kutta) \n');
disp('____________________________');
fprintf(' %5.4f %11.8f\n',x(1),y(1));
for i = 1:6
x(i+1) = x(i)+h;
k1 = f(x(i), y(i));
k2 = f(x(i)+h/12, y(i) + h*k1/12);
k3 = f(x(i)+h/6, y(i)+h/6*(k2));
k4 = f(x(i)+h/4, y(i)+h/16*(k1+3*k3));
k5 = f(x(i)+3*h/4, y(i)+h/16*(21*k1 - 81*k3 + 72*k4));
k6 = f(x(i)+16*h/17, y(i)+h/250563*(1344688*k1...
- 5127552*k3 + 4096896*k4 - 78208*k5));
k7 = f(x(i)+h/2, y(i)+h/234624*( -341549*k1 + 1407744*k3...
- 1018368*k4 + 84224*k5 - 14739*k6));
k8 = f(x(i)+h, y(i)+h/136864*( -381875*k1 + 1642368*k3...
- 1327872*k4 + 72192*k5 + 14739*k6 + 117312*k7));
k9 = f(x(i)+2*h/3, y(i)+h/16755336*( -2070757*k1 + 9929088*k3...
+ 584064*k4 + 3023488*k5 - 447083*k6 + 151424*k7));
k10= f(x(i)+h, y(i)+h/10743824*( 130521209*k1 - 499279872*k3...
- 391267968*k4 + 13012608*k5 - 3522621*k6 + 9033024*k7...
- 30288492*k9));
y(i+1) = y(i)+h/90*(7*k1+32*k4+32*k5+12*k7+7*k8);
fprintf(' %5.4f %11.8f\n',x(i+1),y(i+1));
plot(x,y,'linewidth',2); hold on;
end
for i = 7:n
xx0 = x0+i*h;
pred = y(7) + h*(198721*f(x(7),y(7)) - 447288*f(x(6),y(6))...
+ 705549*f(x(5),y(5)) - 688256*f(x(4),y(4)) +...
407139*f(x(3),y(3))- 134472*f(x(2),y(2)) +...
19087*f(x(1),y(1)))/60480;
pred = y(7) + h*(19087*f(xx0,pred) + 65112*f(x(7),y(7))...
- 46461*f(x(6),y(6)) + 37504*f(x(5),y(5)) - ...
20211*f(x(4),y(4))+ 6312*f(x(3),y(3)) -...
863*f(x(2),y(2)))/60480;
fprintf(' %5.4f %11.8f\n',xx0,pred);
for j=1:6
x(j)=x(j+1);
y(j)=y(j+1);
end
x(7)=xx0;
y(7)=pred;
plot(x,y,'-o','linewidth',2);
end
hold off;
disp(['jumlah iterasi = ',num2str(i)]);
% fprintf('Hasil y dengan Metode Adam-Bashforth orde 7= %11.10f\n'...
% ,pred);
disp(['Hasil y dengan Metode Adam-Bashforth orde 7= ', num2str(pred)]);
xlabel('Independent Variable x'); ylabel('Y Value')
title('Adam-Bashforth-Method'); grid on
This example is my successful program (another program and similar to my work) and shows the output correctly. But it just has one ODE. It's written in another languange (Not English), but i think it's no a big deal bcz i've mentioned about the input. This Example might be useful for you to comparing my work above . Please do a correction. There are several errors in my first script.
And again, thank you.
  댓글 수: 1
Looky
Looky 2019년 12월 13일
Hey Shandy,
have a look at your defintion of
fS = @(t,S,E,I,R)(-(0.5*I + 0.25)*S);
but in the ADM part you use for example:
fS(t(7),S(7))
which is not following your defintion of fS with 5 arguments. You can fix that by adding the other arguments to the call in ADM.

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

답변 (0개)

카테고리

Help Center 및 File Exchange에서 Chassis Systems에 대해 자세히 알아보기

태그

제품


릴리스

R2014a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by