Ch21.1 Numerical Methods for First-Order ODEs

In this section we develop numerical methods to solve first order initial value problems (IVPs).

A first order IVP consists of an ordinary differential equation (ODE) together with an initial condition:

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

An example that we will look at is

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

The numerical methods presented in Ch21.1 are Euler's method, Heun's method, and the Runge-Kutta method.

What do you call a belt made out of watches?

A waist of time.

The method of integrating factors can be used to solve our IVP:

\[ \begin{aligned} y' &= x+y, \,\, y(0) = 0 \\ y' - y &= x, \,\, y(0) = 0 \\ \\ \mu(x) & = e^{\int -1 dx} = e^{-x} \\ y(x) & = \frac{\int x e^{-x} dxt}{e^{-x}} + \frac{C}{e^{-x}}, \,\, y(0) = 0 \\ y(x) & = \frac{-x e^{-x}+\int e^{-x}dx}{e^{-x}} + \frac{C}{e^{-x}}, \,\, y(0) = 0 \\ y(x) & = e^x - x - 1 \end{aligned} \]

On grid points below, plot the direction field using slope formula:

\[ y' = x+y \]

Use initial condition to plot solution in slope field:

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

https://www.geogebra.org/m/EGGQ4bjp

A numerical solution uses ODE slope formula \( y' = f(x,y) \) to generate data points that approximate solution; see weblink.

• Suppose \( g(x) \) is an exact solution to the IVP

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

• To generate a numerical solution, we seek values \( y_n \) that approximate \( g(x_n) \) at nodes \( x_n \).
  
• Three vectors of interest are

\[ \begin{aligned} \mathbf{x} &= [x_0,x_1, \ldots, x_{N}] \\ \mathbf{y} &= [y_0,y_1, \ldots, y_{N}] \\ \mathbf{g} &= [g(x_0),g(x_1), \ldots, g(x_{N})] \\ x_n & = a + nh, \,\,\, n = 0, 1, \ldots, N, \,\,\, h = \frac{b-a}{N} \end{aligned} \]

• At each node there will be some error:
  

\[ e_n = g(x_n) - y_n \]

• The error vector is given by \( \mathbf{E} = [e_1, e_2, \ldots, e_N] \).

• The vector 2-norm \( ||\mathbf{E}|| \) measures overall error:

\[ ||\mathbf{E}|| = \sqrt{e_1^2 + e_2^2 + \cdots + e_N^2} \]

• The percent relative error PRE can also be computed:

\[ PRE = \frac{||\mathbf{E}||}{||\mathbf{g}||}\times 100 \]

Euler's method uses tangent line to obtain \( y_{n+1} \) from \( y_n \).

Euler's method uses tangent line to obtain \( y_{n+1} \) from \( y_n \):

\[ y_{n+1} = y_n + hf(x_n,y_n), \,\, y_0 = y(x_0) \]

If we let \( b_n = hf(x_n,y_n) \), then

\[ \begin{array} {rcl} y_1 & = & y_0 + b_0 \\ y_2 & = & y_1 + b_1 \\ \vdots & = & \vdots \\ % y_{n+1} & = & y_n + b_n \\ % \vdots & = & \vdots \\ y_{N+1} & = & y_N + b_N \end{array} \]

• Starting with given initial point \( (x_0,y_0) \), the tangent line directs us to \( (x_1,y_1) \), where \( y_1 = y_0 + hf_0 \), and so on.
  
• For nonlinear solution, \( (x_1,y_1) \) will not be on solution curve, resulting in an error value \( e_1 = g(x_1) - y_1 \).

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

def euler(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 
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 

#Euler Loop
for n in range(0,N):    
    x[n+1] = x[n] + h              #Next x-value 
    y[n+1] = y[n] + h*f(x[n],y[n]) #Next y value 
    gn[n+1] = g(x[n+1])            #Next gn value
    en[n+1] = g(x[n+1])-y[n+1]     #Next en value

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

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

#Plot the numerical and analytical solutions
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','Euler'),loc = 0)
plt.show()

The following command runs Euler's method for \( N = 5 \).

euler(5)

   Euler Solution  Exact Solution  Error Values
0         0.00000        0.000000      0.000000
1         0.00000        0.021403      0.021403
2         0.04000        0.091825      0.051825
3         0.12800        0.222119      0.094119
4         0.27360        0.425541      0.151941
5         0.48832        0.718282      0.229962
Norm of Error Vector =  0.29659857421909624
Percent Relative Error =  34.1292251762755

• Heun's method is a predictor-corrector method.
• First, Euler's method is used to compute a predicted \( y_{n+1}^* \)

\[ \small{ \begin{aligned} k_1 & = hf(x_n,y_n) \\ y_{n+1}^* &= y_n + k_1 \end{aligned}} \]

• Average of predicted updates produces corrected \( y_{n+1} \)

  \[ \small{ \begin{aligned} k_1 & = hf(x_n,y_n) \\ k_2 & = hf(x_n + h,y_n + k1) \\ y_{n+1} &= y_n + (k_1 + k_2)/2 \end{aligned}} \]

• See graphs on next slide (weblink).

Heun's method is a predictor-corrector method.

\[ \small{ \begin{aligned} k_1 & = hf(x_n,y_n) \\ k_2 & = hf(x_n + h,y_n + k1) \\ y_{n+1} &= y_n + (k_1 + k_2)/2 \end{aligned}} \]

Heun's method:

\[ \small{ \begin{aligned} k_1 & = hf(x_n,y_n) \\ k_2 & = hf(x_n + h,y_n + k1) \\ y_{n+1} &= y_n + (k_1 + k_2)/2 \end{aligned}} \]

for n in range(0,N):    
    x[n+1] = x[n] + h    
    k1 = h*f(x[n], y[n])        
    k2 = h*f(x[n] + h, y[n] + k1)  
    y[n+1] = y[n] + (k1 + k2)/2 

The following command runs Heun's method for \( N = 5 \).

heun(5)

   Heun Solution  Exact Solution  Error Values
0       0.000000        0.000000      0.000000
1       0.020000        0.021403      0.001403
2       0.088400        0.091825      0.003425
3       0.215848        0.222119      0.006271
4       0.415335        0.425541      0.010206
5       0.702708        0.718282      0.015574
Norm of Error Vector =  0.01999320457886825
Percent Relative Error =  2.300592923833429

• RK4 is a predictor-corrector method, using 4 predictors.
  
• First, a full step predictor is computed, followed by two half-step predictors, then another full step predictor, and finally the corrector value for \( y_{n+1} \):

\[ \begin{array} {rcl} k_1 & = & hf(x_n,y_n) \\ k_2 & = & hf(x_n+0.5h,y_n+0.5k_1) \\ k_3 & = & hf(x_n+0.5h,y_n+0.5k_2) \\ k_4 & = & hf(x_n+h,y_n+k_3) \\ y_{n+1} & = & y_n + \frac{1}{6}\left[ k_1+2k_2+2k_3+k_4\right] \end{array} \]

• See next slide for illustration.

\[ \begin{array} {rcl} k_1 & = & hf(x_n,y_n) \\ k_2 & = & hf(x_n+0.5h,y_n+0.5k_1) \\ k_3 & = & hf(x_n+0.5h,y_n+0.5k_2) \\ k_4 & = & hf(x_n+h,y_n+k_3) \\ y_{n+1} & = & y_n + \frac{1}{6}\left[ k_1+2k_2+2k_3+k_4\right] \end{array} \]

for n in range(0,N): 
    x[n+1] = x[n] + h  
    k1 = h*f(x[n],y[n])        
    k2 = h*f(x[n]+0.5*h,y[n]+0.5*k1)  
    k3 = h*f(x[n]+0.5*h,y[n]+0.5*k2)      
    k4 = h*f(x[n+1],y[n]+k3)          
    y[n+1] = y[n] + (k1+2*k2+2*k3+k4)/6  

https://www.geogebra.org/m/b8makdud)

The following command runs RK4 for \( N = 5 \).

rk4(5)

   RK4 Solution  Exact Solution  Error Values
0      0.000000        0.000000      0.000000
1      0.021400        0.021403      0.000003
2      0.091818        0.091825      0.000007
3      0.222106        0.222119      0.000012
4      0.425521        0.425541      0.000020
5      0.718251        0.718282      0.000031
Norm of Error Vector =  3.9388928588386105e-05
Percent Relative Error =  0.004532434509453257

https://www.geogebra.org/m/hvEzUcDF

• For \( N = 5 \) and \( 20 \), the predictor-corrector approach taken by Heun's method is more accurate than Euler's method, while the RK4 method is the most accurate of the three.
• For \( N=5 \), the results, errors and graphs for each numerical method matches those listed in the book (Table 1.1 & Figure 9 on page 11; Tables 21.2, 21.4, 21.5 in Ch21.1).
• The pointwise errors \( e_n \) for each method decrease when using \( N=20 \) nodes compared to using \( N=5 \) nodes.
• Similarly for error norm \( || \mathbf{E} || \) and percent relative error.
• See next slide for comparison of PREs.

• The tables below list PRE ratios for the three methods.
  
• From these ratios, we see that there is a signicant jump in accuracy as the complexity of the method increases.
• The RK4 method is a good blend of accuracy and complexity.

  PRE Ratios  Euler/Heun      Euler/RK4
0      N = 5   14.782609    7555.555556
1     N = 20   58.823529  476190.476190

  PRE Ratios     Heun/RK4
0      N = 5   511.111111
1     N = 20  8095.238095