Two Step Adam Bashford Method
47 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I am trying to make a function that implements the two step Adam Bashford Method to solve an ODE
function [t, w, h] = abs2(f, a, b, alpha, n)
%AB2 Two-step Adams Bashforth method
% [t, w, h] = ab2(f, a, b, alpha, n) performs the two-step Adams Bashforth
% method for solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b. The first step from t = a to t = a + h
% is performed using the modified Euler method.
h=(b-a)/n;
t=a:h:b;
w=zeros(1,length(t));
w(1)=alpha;
%modified euler method
for i=2
k1=h*f(t(i),w(i));
k2=h*f(t(i)+h,w(i)+k1);
w(i+1)=w(i)+1/2*(k1+k2);
end
for i=3:length(t)
w(i+1)=w(i)+(3/2)*h*f(t(i),w(i))-.5*h*f(t(i-1),w(i-1));
end
end
Code to call the function
f=@(t,y) 3*t+y/t;
alpha=5;
a=1;
b=2;
n=3;
[t, w, h] = ab2(f, a, b, alpha, n)
%% Output %%
t = 1×4
1.0000 1.3333 1.6667 2.0000
w = 1×5
5.0000 0 1.6333 3.9567 6.9492
h = 0.3333
As you can the output doesn't look right at all. Any help would be appreciated
0 commentaires
Réponse acceptée
Alan Stevens
le 21 Sep 2021
As follows
f=@(t,y) 3*t+y/t;
alpha=5;
a=1;
b=2;
n=3;
[t, w, h] = abs2(f, a, b, alpha, n);
plot(t,w,'-o'),grid
xlabel('t'),ylabel('w')
disp(['h = ' num2str(h)])
function [t, w, h] = abs2(f, a, b, alpha, n)
%AB2 Two-step Adams Bashforth method
% [t, w, h] = ab2(f, a, b, alpha, n) performs the two-step Adams Bashforth
% method for solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b. The first step from t = a to t = a + h
% is performed using the modified Euler method.
h=(b-a)/n;
t=a:h:b;
w=zeros(1,length(t));
w(1)=alpha;
%modified euler method
k1=h*f(t(1),w(1));
k2=h*f(t(1)+h,w(1)+k1);
w(2)=w(1)+1/2*(k1+k2);
for i=2:length(t)-1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w(i+1)=w(i)+(3/2)*h*f(t(i),w(i))-.5*h*f(t(i-1),w(i-1));
end
end
Plus de réponses (1)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!