How to make convergence plot (error VS time step) in a log-log scale among four numerical methods and exact solution?

조회 수: 41 (최근 30일)
Greetings!
Can someone help and guide me to make a convergence plot (error VS time step) in a log-log scale among four numerical methods and exact solution? The numerical methods are: 1st-order Adams-Bashforth (Explicit), 2nd-order Adams-Bashforth (Explicit), 1st-order Adams-Moulton (Implicit), and 2nd-order Adams-Moulton (Implicit). Here I attached my script, you may review the code. The line should be five, i.e., four are the numerical methods' solutions and one is the exact solution.
%% 1st-order Adams-Bashforh Solution
fun = @(t,y) (y); %Function f(t,y)
y0 = 1,2; %Initial Condition
a = 0; %Starting Time
b = 10; %Ending Time
N = 5; %Number of subintervals, h = (b-a)/N
[t,y] = eul(fun,a,b,y0,N);
% Display Solution
exactY = @(t) exp(t);
Y = exactY(t);
disp('---------------------------------------------------------')
disp([' Euler Method (1st-order Adams-Bashforth); h= ', num2str((b-a)/N)])
disp('t_i y_i y(t_i) |y_i-y(t_i)|')
disp('---------------------------------------------------------')
formatSpec = '%.4f %.4f %.4f %.4f\n';
fprintf(formatSpec,[t';y';Y';(abs(y-Y))'])
%% 2nd-order Adams-Bashforh Solution
fun = @(t,y) (y);
y0 = 1,2;
tspan = [0,10];
N = 4;
% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t) exp(t);
t1 = tspan(1) + h; y1 = exactY(t1);
t2 = tspan(1) + 2*h; y2 = exactY(t2);
[t2,Y2] = AB2(fun,tspan,y0,y1,N);
% Display Solution
Y = exactY(t2);
plot(t2,Y2,t2,Y)
legend('2nd-order Adams-Bashforth Solution','Analytical Solution')
grid on
disp('-----------------------------');
disp('t_i y(t_i) AB2 Error')
disp('-----------------------------');
formatSpec = '%.2f %.5f %.5f %.5f\n';
fprintf(formatSpec,[t2';Y';Y2';abs(Y'-Y2')])
%% 1st-order Adams-Moulton Solution
fun = @(t,y) (y); %Function f(t,y)
y0 = 1,2; %Initial Condition
a = 0; %Starting Time
b = 10; %Ending Time
N = 4; %Number of subintervals, h = (b-a)/N
[t,y] = AM1(fun,a,b,y0,N);
% Display Solution
exactY = @(t) exp(t);
Y = exactY(t);
disp('---------------------------------------------------------')
disp([' 1st-order Adams-Moulton; h= ', num2str((b-a)/N)])
disp('t_i y_i y(t_i) |y_i+y(t_i)|')
disp('---------------------------------------------------------')
formatSpec = '%.4f %.4f %.4f %.4f\n';
fprintf(formatSpec,[t';y';Y';(abs(y-Y))'])
%% 2nd-order Adams-Moulton Solution
fun = @(t,y) (y);
y0 = 1,2;
tspan = [0,10];
N = 4;
% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t) exp(t);
t1 = tspan(1) + h; y1 = exactY(t1);
t2 = tspan(1) + 2*h; y2 = exactY(t2);
[t2,Y2] = AM2(fun,tspan,y0,y1,N);
% Display Solution
disp('-----------------------------');
disp('t_i y(t_i) AM2 Error')
disp('-----------------------------');
formatSpec = '%.2f %.5f %.5f %.5f\n';
fprintf(formatSpec,[t2';Y';Y2';abs(Y'-Y2')])
%% Plot All the Solutions
% Plot Solution of AB1
fplot(exactY, [a,b],'r-','linewidth',2)
hold on
plot(t,y,'b*-','linewidth',2)
hold on
legend('Analytical Solution','Euler (1st-order Adams-Bashforth) Solution','Location','northwest')
% Plot Solution of AB2
fplot(exactY, [a,b],'r-','linewidth',2)
hold on
plot(t,y,'b*-','linewidth',2)
hold on
legend('Analytical Solution','1st-order Adams-Moulton Solution','Location','northwest')
% Plot Solution of AM1
fplot(exactY, [a,b],'r-','linewidth',2)
hold on
plot(t,y,'b*-','linewidth',2)
hold on
legend('Analytical Solution','1st-order Adams-Moulton Solution','Location','northwest')
% Plot Solution of AM2
Y = exactY(t2);
plot(t2,Y2,t2,Y)
legend('2nd-order Adams-Moulton Solution','Analytical Solution')
hold off
grid on
That's all. Please let me know if you need the code of the functions. I really appreciate any help you can provide. Thank you so much.
  댓글 수: 2
Alan Stevens
Alan Stevens 2022년 3월 6일
You'll need something like
errorAM = exacty - Y2;
for each numerical routine, then use the loglog function to do the plotting.
Yusuf Arya Yudanto
Yusuf Arya Yudanto 2022년 3월 7일
Hello. Thank you for your reply. Do you mind to show me the code for the error and loglog function (only for one numerical method)? I really appreciate your help. Thank you
function [t,y] = AM2(fun,tspan,y0,y1,N)
a = tspan(1);
b = tspan(2);
h = (b-a)/N;
t = zeros(N+1,1);
y = zeros(N+1,1);
t(1) = a; y(1) = y0;
t(2) = a+h; y(2) = y1;
for i = 2:N
t(i+1) = t(i) + h;
Fi = fun(t(i),y(i));
Fi1 = fun(t(i-1),y(i-1));
funh = @(yh) yh - y(i) - h/12*(5*fun(t(i+1),yh) + 8*Fi - Fi1);
y(i+1) = fsolve(funh,y(i));
end
%% 2nd-order Adams-Moulton Solution
fun = @(t,y) ((1+4*t)*((y)^1/2));
y0 = 1;
tspan = [0,10];
N = 5;
% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t)(t/2 + t.^2 + 1).^2;
t1 = tspan(1) + h; y1 = exactY(t1);
t2 = tspan(1) + 2*h; y2 = exactY(t2);
[t2,Y2] = AM2(fun,tspan,y0,y1,N);
% Display Solution
disp('-----------------------------');
disp('t_i y(t_i) AM2 Error')
disp('-----------------------------');
formatSpec = '%.2f %.5f %.5f %.5f\n';
fprintf(formatSpec,[t2';Y';Y2';abs(Y'-Y2')])

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

답변 (1개)

Alan Stevens
Alan Stevens 2022년 3월 7일
This should give you the right idea:
%% 1st-order Adams-Bashforh Solution
fun = @(t,y) (y); %Function f(t,y)
y0 = 1; %Initial Condition
a = 0; %Starting Time
b = 10; %Ending Time
N = 5; %Number of subintervals, h = (b-a)/N
[t,y] = eul(fun,a,b,y0,N);
% Display Solution
exactY = @(t) exp(t);
Y = exactY(t);
err = y - Y;
loglog(t,err,'o--')
function [t,yeul] = eul(fun, a, b, y0, N)
h = (b-a)/N;
yeul(1) = y0;
t(1) = 0;
for i = 2:N
t(i) = t(i-1) + h;
yeul(i) = yeul(i-1)+h*fun(t(i-1),yeul(i-1));
end
end

카테고리

Help CenterFile Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by