Ch21.2 Multistep Methods

• Ch21.2 introduces the Adams-Moulton method.
• This method is used to find the numerical solution of the IVP

\[ y' = f(x,y), \,\, y(x_0) = x_0 \]

• As in Ch21.1, the IVP we consider is

\[ y' = x+y, \,\, y(0) = 0 \]

• A one-step method uses the previous \( y_n \) to compute \( y_{n+1} \).
• A multi-step method two or more previous values.

Wanna hear a joke about paper?

Never mind-it's tearable.

• The Adams-Moulton method is a multistep implicit predictor-corrector method that uses four values from previous steps
  • Starting values \( y_1, y_2, y_3 \) obtained using RK4
  • Adams-Bashforth formula used for predicted value \( y_{n+1}^* \)
  • Adams-Moulton formula used for corrected value \( y_{n+1} \)

\[ \small{ \begin{array} {rcl} y' & = & f(x,y), \,\, y(x_0) = y_0 \\ y_{n+1}^* & = & y_n + \frac{h}{24}\left(55f_n -59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\ y_{n+1} & = & y_n + \frac{h}{24}\left(9f_{n+1}^* + 19f_n - 5f_{n-1} + f_{n-2}\right) \end{array} } \]

• Adams-Moulton is a non-self starting and implicit method
  • Use RK4 and \( y_0 \) to obtain \( y_1,y_2,y_3 \) (Red)
  • Use Adams-Bashforth-Moulton after that (Blue)

\[ \small{ \begin{array} {rcl} y_{n+1}^* & = & y_n + \frac{h}{24}\left(55f_n -59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\ y_{n+1} & = & y_n + \frac{h}{24}\left(9f_{n+1}^* + 19f_n - 5f_{n-1} + f_{n-2}\right) \end{array} } \]

• The Adams-Bashforth formula for \( y_{n+1} \) is obtained from using a cubic interpolation polynomial \( p_3(x) \) to fit
  

\[ \small{(x_n,f_n), (x_{n-1},f_{n-1}), (x_{n-2},f_{n-2}), (x_{n-3},f_{n-3})} \]

• Recall here that \( y'=f(x,y) \) and \( f_n = f(x_n,y_n) \).

• The Adams-Moulton formula for \( y_{n+1} \) is obtained from a cubic interpolation polynomial \( p_3(x) \) to fit

\[ \small{(x_{n+1},f_{n+1}), (x_{n},f_{n}), (x_{n-1},f_{n-1}), (x_{n-2},f_{n-2})} \]

• A Newton interpolation method is used in the reading, but a Lagrange approach can be used as well.

\[ \small{ \begin{aligned} p_1(x) &= \frac{x-x_1}{x_0-x_1}f_0 + \frac{x-x_0}{x_1-x_0}f_1 \\ y_{n+1} & = y_n + \int_{x_n}^{x_{n+1}} p_1(x) dx \\ & = y_n + \frac{h}{2} (f_0 + f_1) \,\, (\mathrm{Maple}) \end{aligned} } \]

\[ \small{ y_{n+1}^* = y_n + \frac{h}{24}\left(55f_n -59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) } \]

\[ \small{ y_{n+1} = y_n + \frac{h}{24}\left(9f_{n+1}^* + 19f_n - 5f_{n-1} + f_{n-2}\right) } \]

• Use RK4 and \( y_0 \) to obtain \( y_1,y_2,y_3 \) (Red)
• Use Adams-Bashforth-Moulton formulas after that (Blue)

\[ \small{ \begin{array} {rcl} y_{n+1}^* & = & y_n + \frac{h}{24}\left(55f_n -59f_{n-1} + 37f_{n-2} - 9f_{n-3}\right) \\ y_{n+1} & = & y_n + \frac{h}{24}\left(9f_{n+1}^* + 19f_n - 5f_{n-1} + f_{n-2}\right) \end{array} } \]

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import norm
import pandas as pd

def adamsm(N):  
  #Derivative function from ODE
  f = lambda x,y: x + y    
  #Exact solution; see page 10.
  g = lambda x: np.exp(x)-x-1.0 

a = 0    #Left endpoint of [a,b]
b = 1    #right endpoint of [a,b] 
h = (b-a)/N   #Compute step size h  
x0 = a        #Initial x value as given in IVP
y0 = 0        #Initial y value as given in IVP
x = np.zeros(N+1)     #Initialize x values 
y = np.zeros(N+1)     #Initialize y values 
gn = np.zeros(N+1)    #Initialize gn values
en = np.zeros(N+1)    #Initialize en values
x[0] = x0             #Initial x value 
y[0] = y0             #Initial y value 
gn[0] = g(x0)         #Initial gn value 
en[0] = g(x0) - y0    #Initial en value 

# Will need code here (Homework)
# This includes an initial loop for RK4
#    - RK4 computes the starter values
#    - Updates for gn and en also
# After that, will need a loop for Adams-Moulton
#    - Adams-Bashforth predictor
#    - Adams-Moulton corrector
#    - Updates for gn and en also

#Create a table for iterates    
table = pd.DataFrame({'Adams-Moulton':y,'Exact Solution':gn,'Error Values':en})
print(table) 

NormE = norm(en)    
print('Norm of Error Vector = ', NormE )
PRE = NormE/norm(gn)*100  
print('Percent Relative Error = ', PRE )

L = N*h  #[x0,x0+L] is x-interval for solution
xvalues = np.linspace(x0,x0+L,1000) 
plt.figure(figsize=(12,6))
plt.plot(xvalues,g(xvalues),'r',x,y,'bo') 
plt.xlabel('x')
plt.legend(('Exact','Adams-Moulton'),loc = 0)
plt.show()

For \( N = 5 \), we obtain the following results.

adamsm(5)

   Adams-Moulton  Exact Solution  Error Values
0       0.000000        0.000000      0.000000
1       0.091733        0.091825      0.000091
2       0.425268        0.425541      0.000273
3       1.119507        1.120117      0.000610
4       2.351724        2.353032      0.001309
5       4.386499        4.389056      0.002557
Norm of Error Vector =  0.0029508911369835384
Percent Relative Error =  0.05760125375261964

adamsm(10)

    Adams-Moulton  Exact Solution  Error Values
0        0.000000        0.000000  0.000000e+00
1        0.021400        0.021403  2.758160e-06
2        0.091818        0.091825  6.737641e-06
3        0.222106        0.222119  1.234405e-05
4        0.425528        0.425541  1.305017e-05
5        0.718269        0.718282  1.313731e-05
6        1.120104        1.120117  1.276326e-05
7        1.655188        1.655200  1.156077e-05
8        2.353023        2.353032  9.194714e-06
9        3.249642        3.249647  5.215050e-06
10       4.389057        4.389056 -9.774843e-07
Norm of Error Vector =  3.094385494631146e-05
Percent Relative Error =  0.0004885916788014132

• For \( N=10 \), our results matches those found in Table 21.9.
  
• The pointwise errors \( e_n \) decrease when using \( N=20 \) nodes compared to using \( N=5 \) nodes.
• Similarly for error norm \( || \mathbf{E} || \) and percent relative error.
• The PRE values for Adams-Moulton are larger than RK4 by a factor of about 14; see error ratios below.

  PRE Ratios     AM/RK4
0      N = 5  12.800000
1     N = 20  14.761905