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

12 views (last 30 days)
Yusuf Arya Yudanto
Yusuf Arya Yudanto on 5 Mar 2022
Answered: Alan Stevens on 7 Mar 2022
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
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

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!