How can I implement Richardson Extrapolation to the third order Adams-Bashforth Method

Question

Anil Aydin Aytekin 2020-11-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/669738-how-can-i-implement-richardson-extrapolation-to-the-third-order-adams-bashforth-method

在 MATLAB Online 中打开

So, there is an second order differantial equation of ;

which has the initial condition of;

. And the analytical solution is;

.

the desired local accuracy equal to 10^-4

where the accuracy is defined with:

So we need to expect the condition below;

so far, I have done the analytical solution and Adams-Bashforth 3rd method but could not understand how to implement the richardson method to the Adams method. And I need to get a result of a single plot with two traces: analytical solution and numerical solution satisfying the conditions above.

I know the mathematical formulas of richardson method, but i am not familiar it into matlab. If anyone can help me on this problem, that would be really great.

Thanks all.

These are codes that I have wrote ;

%%% Analytical solution %%%%%
h  = 0.1;
x  = 0:h:5;
c1 = 1;
c2 = (1 - (10*pi)/(1000*(pi)^2 -1 )) / sqrt(10);
y_analytical = c2*sin( sqrt(10).*x ) + c1*cos( sqrt(10).*x )+ sin( 100*pi.*x )/( 10000*pi^2-10 );    
plot (x,y_analytical);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  The code below is a seperate script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ADAM'S BASFORTH METHOD %%%
clear all
h = 0.01;
y = [1;1];
t = 0:h:0.2;
y_analitical = funct_analitical(t);
y(:,2) = y(:,1)+h*funct(t(1),y(:,1));
y(:,3) = y(:,2) + 1/2*h*(3*funct(t(2), y(:,2)) - funct(t(2)-h, y(:,1)));
for i = 3 : length(t)-1
    y_minus2 = [
        spline(t(1:i),y(1,:), t(i)-2*h)
        spline(t(1:i),y(2,:), t(i)-2*h)
    ];
    y_minus1 = [
        spline(t(1:i),y(1,:), t(i)-h)
        spline(t(1:i),y(2,:), t(i)-h)
    ];
    y(:,i+1) = y(:,i) + h*(23/12*funct(t(i), y(:,i)) - 4/3*funct(t(i)-h, y_minus1) +5/12*funct(t(i)-2*h, y_minus2)); % Second order
end
%%% ERROR CALCULATIONS %%% 
avg_error = mean(abs(y(1,:)- y_analitical))
no_of_steps = length(t)
max_error = max (abs(y(1,:)- y_analitical))
figure
plot(t,y_analitical,'black');
hold on
plot(t,y(1,:),'red');
grid on
%%% ANALİTİCAL SOLUTION %%%
function out = funct_analitical(t)
c1 = 1;
c2 = (1 - (10*pi)/(1000*(pi)^2 -1 )) / sqrt(10);
out = c2*sin( sqrt(10).*t ) + c1*cos( sqrt(10).*t ) + sin( 100*pi.*t )/( 10000*pi^2-10 );
end
%%% EQUATİON %%%
function out = funct(t,y)
dy2dt = -10*y(1,1)-sin(100*pi*t);
dy1dt = y(2,1);
out = [dy1dt; dy2dt];
end