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.
My code below:
The $1e-10$ is quite arbitary. Here is some output:
This matches with the results you require. Note, that I perform one fixpoint iteration more per each step. Maybe my tolerance is too small.