Three step Adams-Moulton functional iteration

Question

I'm given an IVP:

$$y' = e^y, \;\;\;\; 0 \le t \le 0.20, \;\;\;\;\; y(0) = 1$$ The solution is: $$y(t) = 1 - \ln(1 - et)$$

Applying the three-step Adams-Moulton method to the problem is equivalent to finding the fixed point $w_{i + 1}$ of

$$g(w) = w_i + \frac{h}{24}(9e^w + 19e^{w_i} - 5e^{w_{i - 1}} + e^{w_{i - 2}})$$

The directions are to use $h = 0.01$ to obtain $w_{i + 1}$ by functional iteration for $i = 2, ..., 19$ using exact starting values $w_0$, $w_1$, and $w_2$. At each step, use $w_i$ to initially approximate $w_{i + 1}$.

Here's the method as described in the book:

Adams-Moulton Three-Step Implicit Method $$w_0 = \alpha, \;\;\;\; w_1 = \alpha_1, \;\;\;\; w_2 = \alpha_2$$ $$w_{i + 1} = w_i + \frac{h}{24}[9f(t_{i + 1}, w_{i+1}) + 19f(t_i, w_i) - 5f(t_{i-1}, w_{i-1}) + f(t_{i-2}, w_{i - 2})]$$ where $i = 2, 3, ..., N - 1$.

I tried implementing this in MATLAB but I wasn't getting the answers as shown on page 4 of this document. (The book also had the same answers.)

Here's my code:

function adamsMoulton3(a, b, alpha, h, f)
    N = (b - a)/h;
    t = @(i) a + i*h;

    w = zeros(1, N); %Initialize w
    w(1, 1:3) = alpha; %Set initial values w(1), w(2), w(3) to y(0), y(0.01), y(0.02)
    for i = 3:N - 1
        w(i + 1) = w(i) + (h/24)*(9*f(t(i + 1), w(i + 1))...
                        + 19*f(t(i), w(i)) - 5*f(t(i - 1), w(i - 1))...
                        + f(t(i - 2), w(i - 2)));
        fprintf('t = %.2f; w = %.10f\n', t(i + 1), w(i));
    end
end

My commands and output:

>> f = @(t, y) exp(y);
>> y = @(t) 1 - log(1 - exp(1)*t);
>> a = 0; b = 0.20; alpha = [y(0) y(0.01) y(0.02)]; h = 0.01;
>> adamsMoulton3(a, b, alpha, h, f);
t = 0.04; w = 1.0558992926
t = 0.05; w = 1.0777175088
t = 0.06; w = 1.0999020071
t = 0.07; w = 1.1225096281
t = 0.08; w = 1.1455501124
t = 0.09; w = 1.1690419183
t = 0.10; w = 1.1930048097
t = 0.11; w = 1.2174598613
t = 0.12; w = 1.2424295869
t = 0.13; w = 1.2679380734
t = 0.14; w = 1.2940111310
t = 0.15; w = 1.3206764603
t = 0.16; w = 1.3479638423
t = 0.17; w = 1.3759053504
t = 0.18; w = 1.4045355925
t = 0.19; w = 1.4338919837
t = 0.20; w = 1.4640150581

I also tried using the $g(w)$ formula, placing this inside of the for loop:

W = w(i) + (h/24)*(9*exp(W) + 19*exp(w(i)) - 5*exp(i - 1) + exp(w(i - 2)));

I wasn't sure what to use as the first value of W, so I used 1. But that got me vastly different answers.

Is there something wrong with the code?

Answer 1

There are two things wrong in your code.

1. You have to shift $t$ by minus $h$, as for $i=1$ you need to have $t=a$
2. Way more important: For each $w_{i+1}$ you need to perform a fixpoint iteration! This means that you want to have $w_{i+1}$ such that $w_{i+1}=g(w_{i+1})$. That's why the function is defined as $g(w)$.

My code below:

function adamsMoulton3(a, b, alpha, h, f)
    N = (b - a)/h;
    t = @(i) a + (i-1)*h;
    g = @(wFixPoint,w,f,i) w(i) + (h/24)*(9*f(t(i + 1), wFixPoint)...
                    + 19*f(t(i), w(i)) - 5*f(t(i - 1), w(i - 1))...
                    + f(t(i - 2), w(i - 2)));
    % note, that i define g here, as this makes the following 
    % more readable


    alpha
    w = zeros(1, N); %Initialize w
    w(1, 1:3) = alpha; %Set initial values w(1), w(2), w(3) to y(0), y(0.01), y(0.02)
    for i = 3:N-1
        % we need to perform a fixpoint iteration for w(i+1)
        % to get w_i+1 = g(w_i+1)
        % start fixpoint iteration for w_i+1 with w_i as described in a)
        wFixPoint = w(i);
        while (abs(wFixPoint-g(wFixPoint,w,f,i))> 1e-10)
                wFixPoint = g(wFixPoint,w,f,i);
                fprintf('\t doing a fixpointiteration... w =%.10f\n',wFixPoint)
        end
        w(i+1) = wFixPoint;    
        fprintf('t = %.2f; w = %.10f\n', t(i+1), w(i+1));
    end
end

The $1e-10$ is quite arbitary. Here is some output:

         doing a fixpointiteration... w =1.0847471052
         doing a fixpointiteration... w =1.0850626018
         doing a fixpointiteration... w =1.0850661028
         doing a fixpointiteration... w =1.0850661417
         doing a fixpointiteration... w =1.0850661421
t = 0.03; w = 1.0850661421
         doing a fixpointiteration... w =1.1147708245
         doing a fixpointiteration... w =1.1151054513
         doing a fixpointiteration... w =1.1151092778
         doing a fixpointiteration... w =1.1151093215
         doing a fixpointiteration... w =1.1151093220
t = 0.04; w = 1.1151093220
         doing a fixpointiteration... w =1.1457233327
         doing a fixpointiteration... w =1.1460788838
         doing a fixpointiteration... w =1.1460830774
         doing a fixpointiteration... w =1.1460831269
         doing a fixpointiteration... w =1.1460831275
t = 0.05; w = 1.1460831275
         doing a fixpointiteration... w =1.1776638941
         doing a fixpointiteration... w =1.1780423953
         doing a fixpointiteration... w =1.1780470046
         doing a fixpointiteration... w =1.1780470607
         doing a fixpointiteration... w =1.1780470614
t = 0.06; w = 1.1780470614
         doing a fixpointiteration... w =1.2106576289
         doing a fixpointiteration... w =1.2110613761
         doing a fixpointiteration... w =1.2110664578
         doing a fixpointiteration... w =1.2110665218
         doing a fixpointiteration... w =1.2110665226
t = 0.07; w = 1.2110665226
         doing a fixpointiteration... w =1.2447763102
         doing a fixpointiteration... w =1.2452079161
         doing a fixpointiteration... w =1.2452135371
         doing a fixpointiteration... w =1.2452136103
         doing a fixpointiteration... w =1.2452136112
t = 0.08; w = 1.2452136112
...
...
t = 0.19; w = 1.7266507925 % the last result

This matches with the results you require. Note, that I perform one fixpoint iteration more per each step. Maybe my tolerance is too small.