MATLAB Error Implementing Fourth Order Adams-Moulton Method

Question

Trying to implement the fourth order AM method in MATLAB using fourth order RK to get the first four starting values. Code is as follows:

% Use the fourth order Adams-Moulton (AM) method with h = 0.1

% Define the exact solution
y = @(t) (7*exp(3*t)+6*t+2)/9;

% Define the RHS of the differential equation
f = @(t,y) 3*y-2*t;

% Define the step size and create an interval of the solution from 0 to 1
h = 0.1;
t = 0:h:1;

% Initialize the approximate solution, ya
ya = zeros(size(t));

% Use fourth order Runge-Kutte to obtain the first four starting values
for j = 1:4
    ya(1) = f(t(j), y(j));
    ya(2) = f(t(j) + h/2, y(j) + h/2*ya(1));
    ya(3) = f(t(j) + h/2, y(j) + h/2*ya(2));
    ya(4) = f(t(j) + h, y(j) + h*ya(3));
end

% Implement the four step AM method for the remaining time intervals
for j = 1:(length(t)-4)
    ya(j+4) = ya(j+3) + (h/24)*(9*f(t(j+4),y(j+4))+19*f(t(j+3),ya(j+3))-5*f(t(j+2),ya(j+2))+f(t(j+1),ya(j+1)));
end

% Plot the exact solution over the interval [0,1] along with the numerical
% solution for comparison
t1 = linspace(0,1);
plot(t1,y(t1),t,ya,'x','MarkerSize',12,'LineWidth',2)
xlabel('\leftarrow t \rightarrow');ylabel('\leftarrow y \rightarrow');
legend('y(t)-exact solution','y_i - Approx solution, 0\leq i \leq n','Location','NorthWest')
title('Exact solution vs 4^{th} order Adams-Moulton (AM) method')

% Print the numerical solution 
var = {'t','Numerical_Values','Exact_Values'};
T = table(t',ya',y(t)','VariableNames',var);
disp('The result of the fourth order Adams-Moulton (AM) method with step size h=0.1 is:')
disp(T)

Now I'm getting output, but the values of y(1) through y(4) from the RK4 cannot be correct because the values are as follows:

The result of the fourth order Adams-Moulton (AM) method with step size h=0.1 is:

t	Numerical_Values	Exact_Values
0	379769.246644342	1
0.1	436734.533640994	1.33877907255911
0.2	445279.326690492	1.7727590669704
0.3	513352.84465149	2.33524686423318

Thoughts?

Answer 1

You have defined $y$ as an anonymous function

% Define the exact solution
y = @(t) (-2*(-3*t.*exp(-3*t)-exp(-3*t))+7)./(9*exp(-3*t));

This means that you can treat it as a vector

y(j+1) = y(j) + h/6*(k1 + 2*k2 + 2*k3 + k4);

The only way to fix this is to have another vector to store your values per iteration.

Answer 2

On the program structure

The overall problem of this code is that it mixes too many variables with various scopes and functions in one large block of spaghetti code. Mainly it is using y as function and array name, alternatively or in addition using ya as name of the array containing (some of) the computed values.

• You declare the approximate solution as ya as per comment
• You declare y as a function per anonymous function definition
• You try to fill a non-existing array y in the RK4 steps (original code of 2018)
• You proceed to fill the array ya with a zero starter segment

To avoid this general pattern, modularization with procedures and local variables was introduces in times beyond remembering. Now Matlab makes this kind of spaghetti script rather easy, but it also makes an intermediate form of modularization easy, where you structure

main program

function f1

...

function fn

When Matlab reads the script, in a first parsing pass the functions are processed and are available when executing the main program in a second pass.

Main program

So in this case the block of the main program would consist of setting up the ODE and its reference solution, calling the integrator and evaluating its output visually and in its values.

% Define the exact solution
y_exact = @(t) (7*exp(3*t)+6*t+2)/9;

% Define the RHS of the differential equation
f = @(t,y) 3*y-2*t;

% Define the initial point
y0 = y_exact(0)

% Define the step size and create an interval of the solution from 0 to 1
h = 0.1;
t = 0:h:1;

% Integrate 
y_num = AM4int(f,t,y0)

% Plot the exact solution over the interval [0,1] along with the numerical
% solution for comparison
t1 = linspace(0,1,200);
plot(t1,y_exact(t1),t,y_num,'x','MarkerSize',12,'LineWidth',2)
xlabel('\leftarrow t \rightarrow');ylabel('\leftarrow y \rightarrow');
legend('y(t)-exact solution','y_i - Approx solution, 0\leq i \leq n','Location','NorthWest')
title('Exact solution vs 4^{th} order Adams-Moulton (AM) method')

Note that the initial value was taken from the exact solution for consistency, it is easy to check that $y(0)=1$. This information was and is nowhere incorporated in your code.

The Adams-Moulton method

You want to implement an implicit method. Usually the solution of the implicit method step equation is obtained iteratively, by some kind of fixed-point iteration. The classical way is the predictor-corrector method using as predictor the Adams-Bashford method of the same order (see e.g. tables here).

The predictor step is missing in all your code variants, originally there was a variable P(redictor), but without being defined. There is no intention of running a loop for the corrector, so the scheme is a simple PECE (E = evaluation of the ODE function), which makes an explicit method out of the generally implicit one. A fully implicit implementation would run the corrector step until convergence or divergence is detected.

function y = AM4int(f,t,y0)
  h = t(2)-t(1); % equal spacing is assumed in t
  % bootstrap with RK4
  y(1) = y0;
  V(1) = f(t(1),y(1));
  for i = 1:3
    y(i+1) = RK4step(f,t(i),y(i),h);
    V(i+1) = f(t(i+1),y(i+1));
  end%for 

  for i = 4:(len(t)-1)
    % use Adams-Bashford 4 as predictor    
    y(i+1) = y(i) + (h/24)*([-9 37 -59 55]*V');
    % shift the derivative value array, evaluate for the new point
    V(1:3) = V(2:4); V(4) = f(t(i+1),y(i+1));
    % apply Adams-Moulton 4 corrector once  
    y(i+1) = y(i) + (h/24)*([1 -5 19 9]*V');
    % evaluate the derivative for the corrected point
    V(4) = f(t(i+1),y(i+1))
  end%for
end%function

Use the standard RK4 step implementation

function y_next = RK4step(f,t,y,h)
  k1 = h*f(t,y);
  k2 = h*f(t+0.5*h, y+0.5*k1);
  k3 = h*f(t+0.5*h, y+0.5*k2);
  k4 = h*f(t+h, y+k3);
  y_next = y + (k1+2*(k2+k3)+k4)/6;
end%function

The result

t	Numerical_Values	Exact_Values
0	1.0000	1.0000
0.1000	1.3388	1.3388
0.2000	1.7727	1.7728
0.3000	2.3352	2.3352
0.4000	3.0711	3.0712
0.5000	4.0411	4.0413
0.6000	5.3272	5.3275
0.7000	7.0398	7.0404
0.8000	9.3284	9.3291
0.9000	12.3942	12.3953
1.0000	16.5093	16.5110

So indeed the numerical values are correct to the 4 leading digits, giving an exact match visually, as one could expect for a 4th order method with step size $0.1$