How to make graph that plot between Analytical Solution (Exact Solution) and Numerical Method Solution?

Greetings!

Can someone help to guide me how to make graph that plot between Analytical Solution (Exact Solution) and Adams-Moulton Solution? Here I attached my script.

%% 2nd-order Adams-Moulton Solution
fun = @(t,y) ((1+4*t)*((y)^1/2));
y0 = 1;
tspan = [0,1];
N = 4; 
%% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t) (((3*t)/2)+3*t.^2).^(2/3);
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
Y = exactY(t2);
disp('-----------------------------'); 
disp('t_i   y(t_i)   AM2   Error')
disp('-----------------------------'); 
formatSpec = '%.2f   %.5f   %.5f   %.5f\n';
fprintf(formatSpec,[t2';Y';Y2';abs(Y'-Y2')])

I really appreciate any help you can provide. Thank you.

[Torsten]

The Adams-Moulton method is an implicit method.

You can see this by the fact that in your code, Fiplus1 depends on the unknown y(i+1).

As a consequence, you will have to iterate or use a nonlinear equations solver like "fsolve" or "fzero" to solve the equation

y(i+1) = y(i) + (h/12)*(5*fun(t(i+1),y(i+1)) + 8*Fi - Fi1)

for y(i+1).

In your code, you use y(i+1) = 0 on the right-hand side of the equation

y(i+1) = y(i) + (h/12)*(5*fun(t(i+1),y(i+1)) + 8*Fi - Fi1)

because you initialized the y-vector to zero. This is wrong.

[Yusuf Arya Yudanto]

Hi. Thank you for reviewing my script. Do you mind to guide me how to make the right script of Adams-Moulton? I get confused. Where do I supposed to put fsolve or fzero in my script?

[Torsten]

fun = @(t,y) (1+4*t)*sqrt(y);
y0 = 1;
tspan = [0,1];
N = 4; 
%% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t)(t/2 + t.^2 + 1).^2;
t1 = tspan(1) + h; y1 = exactY(t1);
[T,Ynum] = AM2(fun,tspan,y0,y1,N);
%% Display Solution
Yexact = exactY(T);
plot(T,[Ynum,Yexact])

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
end

[Yusuf Arya Yudanto]

Many thanks! I really appreciate your help.

Answer 1

[Alan Stevens]

After Y = exactY(t2) you probably need something like:

plot(t2,Y2,t2,Y)
legend('AM','Exact')

However, you haven't supplied the code for your AM function, so I can't check!

[Yusuf Arya Yudanto]

Well, this is my AM2 function. I'll accept your answer if it succeed.

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));
    Fiplus1 = fun(t(i+1),y(i+1));
    y(i+1) = y(i) + (h/12)*(5*Fiplus1 + 8*Fi - Fi1)
end

Thank you.

[Alan Stevens]

Like so:

%% 2nd-order Adams-Moulton Solution
fun = @(t,y) ((1+4*t)*((y)^1/2));
y0 = 1;
tspan = [0,1];
N = 4; 
%% Initial Values
h = (tspan(2) - tspan(1))/N;
exactY = @(t) (((3*t)/2)+3*t.^2).^(2/3);
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
Y = exactY(t2);

plot(t2,Y2,t2,Y),grid
xlabel('t2'),ylabel('y2 and Y')
legend('AM','Exact')

disp('-----------------------------'); 

-----------------------------

disp('t_i   y(t_i)   AM2   Error')

t_i   y(t_i)   AM2   Error

disp('-----------------------------'); 

-----------------------------

formatSpec = '%.2f   %.5f   %.5f   %.5f\n';
fprintf(formatSpec,[t2';Y';Y2';abs(Y'-Y2')])

0.00   0.00000   1.00000   1.00000
0.25   0.68142   0.68142   0.00000
0.50   1.31037   0.78457   0.52580
0.75   1.99248   0.96652   1.02596
1.00   2.72568   1.26418   1.46150

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));
    Fiplus1 = fun(t(i+1),y(i+1));
    y(i+1) = y(i) + (h/12)*(5*Fiplus1 + 8*Fi - Fi1);
end
end

[Yusuf Arya Yudanto]

Thanks, I got it!

[Alan Stevens]

Looks like you have the exact solution wrong if your function is correct!

Assuming dy/dt = ((1+4*t)*((y)^1/2)) then y = (t/2 + t^2 + 1).^2 (with y(0) = 1)

[Yusuf Arya Yudanto]

Oh. Thank you for the correction. I really appreciate it!

[Yusuf Arya Yudanto]

Hello, Sir Alan. Can you help me make a loglog scale plot of the multiple solutions? I've posted the question, please kindly check it. How to make convergence plot (error VS time step) in a log-log scale among four numerical methods and exact solution?