Adams-Bashforth 4th order

14 visualizaciones (últimos 30 días)
LUCA SCARPELLI
LUCA SCARPELLI el 17 de Jul. de 2021
Editada: LUCA SCARPELLI el 17 de Jul. de 2021
Hi everyone,
I'm trying to implement the AB4 to solve ODE's, and my code seems to be wrong as it's not compiling: I get the error message "Index exceeds the number of array elements (4)."
I tried to check the range of implementation for the code in the AB part (the second loop in the code). I used R-K4 to find the initial condition, as they are four.
Here is my code:
function adams_bashforth_method(a,b,s,f,y0)
% s is the number of subinterval for the interpolation
k=3; % steps-1
x=linspace(a,b, s+1); %network set creation
h=1/(s+1); %discretization of the network set
u=zeros(1,s+1); %network function creation
u(1)=y0; %initial condition of Cauchy problem
for i=1:k
k_1=f(x(i),u(i)); %k_1
k_2=f(x(i)+h./2,u(i)+((h./2).*k_1)); %k_2
k_3=f(x(i)+h./2,u(i)+((h./2).*k_2)); %k_3
k_4=f(x(i)+h,u(i)+(h.*k_3)); %k_4
u(i+1)=u(i)+(h/6).*(k_1+2.*(k_2+k_3)+k_4); %Runge-Kutta
end
%Now that I know the initial four condition, as I set the step constant, I
%can evaluate the others u's by AB4 formula
for i=4:s
u(i+1)=u(i)+(h./24).*(55.*(f(x(i+3),u(i+3)))-59.*(f(x(i+2),u(i+2)))+37.*(f(x(i+1),u(i+1)))-9.*(f(x(i),u(i))));
end
plot(x,u,'r');
title('Adams Bashforth 4th order method')
hold on
end
As a supply, I'll leave here the Cauchy problem I'm facing:
  4 comentarios
Torsten
Torsten el 17 de Jul. de 2021
To dimension x,h and u, you must work with s instead of k.
Further, you must include k_1 = f(x(i),u(i)) in the for-loop in the RK section.
Further, your update in the AB loop for u is wrong:
u(i+1) = u(i) + h/24*(55*f(x(i),u(i))-59*f(x(i-1),u(i-1))+37*f(x(i-2),u(i-2))-9*f(x(i-3),u(i-3)))
Further, you should return u to the calling program from the adams bashforth function.
LUCA SCARPELLI
LUCA SCARPELLI el 17 de Jul. de 2021
ok, that was very useful. Thank you very much!!
I didn't really wrote the return value, as I simply write
function [u]=adams_bashforth_method(a,b,s,f,y0)
It should be the same in this way, isn't it!?

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Productos


Versión

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by