98 ビュー (過去 30 日間)
B 2021 年 9 月 21 日
コメント済み: B 2021 年 9 月 22 日
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)
% [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

サインインしてコメントする。

採用された回答

Alan Stevens 2021 年 9 月 21 日
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)])
h = 0.33333
function [t, w, h] = abs2(f, a, b, alpha, n)
% [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
1 件のコメント表示非表示 なし
B 2021 年 9 月 22 日
thanks

サインインしてコメントする。

カテゴリ

Find more on Programming in Help Center and File Exchange

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by