[Math] Three step Adams-Moulton functional iteration

MATLABnumerical methodsordinary differential equations

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?

Best Answer

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.

Related Question