Computational Physics with Python: Ordinary Differential Equations

Predictor-Corrector and Runge Kutta Methods

Predictor-Corrector Methods

Given any time and state value, the function, F(t,S(t)), returns the change of state dS(t)/dt. Predictor-corrector methods of solving initial value problems improve the approximation accuracy of non-predictor-corrector methods by querying the F function several times at different locations (predictions), and then using a weighted average of the results (corrections) to update the state. Essentially, it uses two formulas: the predictor and corrector. The predictor is an explicit formula and first estimates the solution at t_j+1, i.e. we can use the Euler method or some other methods to finish this step. After we obtain the solution S(t_j+1), we can apply the corrector to improve the accuracy. Using the found S(t_j+1) on the right-hand side of an otherwise implicit formula, the corrector can calculate a new, more accurate solution.

The midpoint method has a predictor step:

which is the prediction of the solution value halfway between tj and tj+1

.

It then computes the corrector step:

which computes the solution at S(t_j+1) from S(t_j) but using the derivative from S(t_j+h2).

Example:

Given the IVP:

import numpy as np
import math 
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt # side-stepping mpl backend
import matplotlib.gridspec as gridspec # subplots
import warnings

warnings.filterwarnings("ignore")

Defining the Function:

def myfun_ty(t,y):
    return t-y

# Start and end of interval
b=2
a=0
# Step size
N=4
h=(b-a)/(N)
t=np.arange(a,b+h,h)
fig = plt.figure(figsize=(10,4))
plt.plot(t,0*t,'o:',color='red')
plt.xlim((0,2))
plt.title('Illustration of discrete time points for h=%s'%(h))

Exact Solution

The initial value problem has the exact solution

The figure below plots the exact solution.

IC=1 # Intial condtion
y=(IC+1)*np.exp(-t)+t-1
fig = plt.figure(figsize=(6,4))
plt.plot(t,y,'o-',color='black')
plt.title('Exact Solution ')
plt.xlabel('time')

IC=1 # Intial condtion
y=(IC+1)*np.exp(-t)+t-1
fig = plt.figure(figsize=(6,4))
plt.plot(t,y,'o-',color='black')
plt.title('Exact Solution ')
plt.xlabel('time')

### Initial conditions
w=np.zeros(len(t))
w0=np.zeros(len(t))
w[0]=IC
w[1]=y[1]
w0[0]=IC
w0[1]=y[1]

for k in range (1,N):
    w0[k+1]=(w[k]+h/2.0*(3*myfun_ty(t[k],w[k])-myfun_ty(t[k-1],w[k-1])))   
    w[k+1]=(w[k]+h/2.0*(myfun_ty(t[k+1],w0[k+1])+myfun_ty(t[k],w[k])))   

Plotting the Solution

def plotting(t,w,w0,y):
    fig = plt.figure(figsize=(10,4))
    plt.plot(t,y, 'o-',color='black',label='Exact')
    plt.plot(t,w0,'v:',color='blue',label='Adams-Bashforth Predictor')
    plt.plot(t,w,'^:',color='red',label='Adams-Moulton Corrector')
    plt.xlabel('time')
    plt.legend()
    plt.show 
plotting(t,w,w0,y)

d = {'time t_i': t, 'Adams Predictor w0': w0,
     'Adams Corrector':w,'Exact (y)':y,'|w-y|':np.round(np.abs(y-w),5),'|w0-w|':np.round(np.abs(w0-w),5),'estimate LTE':np.round(1/(6*h)*abs(w0-w),5)}
df = pd.DataFrame(data=d)
df

Runge Kutta Methods

Runge Kutta (RK) methods are one of the most widely used methods for solving ODEs. Recall that the Euler method uses the first two terms in the Taylor series to approximate the numerical integration, which is linear

We can greatly improve the accuracy of numerical integration if we keep more terms of the series in

In order to get this more accurate solution, we need to derive the expressions of S′′(t_j),S′′′(t_j),⋯,S^(n)(t_j). This extra work can be avoided using the RK methods, which is based on truncated Taylor series, but not require computation of these higher derivatives.

Fourth-order Runge Kutta method

A classical method for integrating ODEs with a high order of accuracy is the Fourth Order Runge Kutta (RK4) method. It is obtained from the Taylor series using the similar approach we just discussed in the second-order method. This method uses four points k1,k2,k3, and k4. A weighted average of these is used to produce the approximation of the solution. The formula is as follows.

Therefore, we will have:

Example:

Solving RK4 for 5x² — y / e^(x+y):

"""Four_Step_Runge_Kutta_ODE1.py 

Implementation of the classic fourth-order method also refered as the
"original" Runge–Kutta method. This method is an implicit four step
Runge-Kutta method which solves an intial value problem numerically. 
"""

from datetime import datetime
import matplotlib.pyplot as plt
from math import exp, sqrt 




def runge_kutta(f, x_0, y_0, h):
    """Four step Runge-Kutta method (RK4)
    Solves first order ODEs
    """
    k_0 = f(x_0, y_0)
    k_1 = f(x_0 + h/2, y_0 + h/2 * k_0)
    k_2 = f(x_0 + h/2, y_0 + h/2 * k_1)
    k_3 = f(x_0 + h, y_0 + h * k_2)

    k = 1/6 * (k_0 + 2.0*k_1 + 2.0*k_2 + k_3)

    x_1 = x_0 + h
    y_1 = y_0 + h * k

    return x_1, y_1


def f(x, y):
    """Example first order ordinary differential equation (ODE)"""
    return (5*x**2 - y) / (exp(x+y))


if __name__=="__main__":
    # Initial values
    x_0 = 0.0
    y_0 = 1.0

    # Step length 
    h = 0.1

    x_values = [x_0]
    y_values = [y_0]

    # Calculate solution
    x = x_0
    y = y_0
    for _ in range(100):
        x, y = runge_kutta(f, x, y, h)
        x_values.append(x)
        y_values.append(y)
        print(x, y)

    # Plot solution
    plt.plot(x_values, y_values)
    plt.show()