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

42 views (last 30 days)
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 Comments
Alan Stevens
Alan Stevens on 6 Mar 2022
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 on 7 Mar 2022
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')])

Sign in to comment.

Answers (1)

Alan Stevens
Alan Stevens on 7 Mar 2022
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

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!