Contents

Tutorial 10 - Initial value problems

Tutorial 10 - Initial value problems

nitial value problems of ordinary differential equations, explicit and implicit Euler method, Runge-Kutta methods, Leap-frog, Adams-Bashford, Adams-Moulton, predictor-corrector, Bulirsch-Stoer algorithm, stiff equations.

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

Runge-Kutta methods

The Runge–Kutta methods are a family of explicit and implicit iterative methods for the numerical solution of the initial value problem

\[ \frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0. \]

Here \( y \) is an unknown function (scalar or vector) of \( x \). The function \( f \) and the initial conditions \( x_0 \), \( y_0 \) are given. The Runge–Kutta methods take the form

\[ y_{i + 1} = y_{i} + h \sum_{j = 1}^{r} b_{j} k_{j}(h), \quad i = 0, \dots, n - 1. \]

where

\[ k_j = f(x_i + c_j h, y_i + h \sum_{l = 1}^{j - 1} a_{jl} k_l) \]

in the case of explicit methods and

\[ k_j = f(x_i + c_j h, y_i + h \sum_{l = 1}^{r} a_{jl} k_l) \]

in the case of implicit methods. Here, \( y_i \) and \( y_{i + 1} \) are the approximations of \( y(x_i)\) and \( y(x_{i + 1}) \), respectively, and \( h > 0 \) is a step size. To specify a particular method, one needs to provide the integer \( r \) (the number of stages), and the coefficients \( a_{ij} \), \( b_i \), and \( c_i \) (see the list of Runge-Kutta methods). The matrix \( (a_{ij}) \) is called the Runge–Kutta matrix, while the \( b_i \) and \( c_i \) are known as the weights and the nodes. These data can be arranged in a Butcher tableau:

\[\begin{split} \begin{array} {c|cccc} c_1 & a_{11} & a_{12} & \cdots & a_{1r} \\ c_2 & a_{21} & a_{22} & \cdots & a_{2r} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_r & a_{r1} & a_{r2} & \cdots & a_{rr}\\ \hline & b_1 & b_2 & \cdots & b_{r} \end{array} \end{split}\]

Note: In the case of explicit methods the Runge-Kutta matrix is lower triangular.

Forward (explicit) Euler method

The forward Euler method

Exercise 10.1: Implement the forward Euler method.

def forward_euler(f, x_0, x_n, y_0, n): # Runge-Kutta 1st order method

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(forward_euler(f, 0, 10, 2, 10000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.2: The fundamental law of radioactive decay can be mathematically expressed as

\[ \frac{dN(t)}{dt} = - \lambda N(t), \qquad N(t_0) = N_0, \]

where \( N \) is the number of radioactive nuclei, \( t \) is time, \( \lambda = \ln(2) \, / \, t_{1/2} \) is the decay constant and \( t_{1/2} \) is the half-life of the decaying quantity. The analytical solution of the decay equation is

\[ N(t) = N_0 e^{-\lambda t}. \]

Solve the decay equation for Ra-226 using the forward Euler method and compare the numerical and analytical solutions.

# add your code here

Runge-Kutta \(2^{\mathrm{nd}}\) order method

Exercise 10.3: Implement the Runge-Kutta \( 2^{\mathrm{nd}} \) order method.

def runge_kutta_2(f, x_0, x_n, y_0, n):

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(runge_kutta_2(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.4: Solve the decay equation for Ra-226 using the Runge-Kutta \( 2^{\mathrm{nd}} \) order method and compare the numerical and analytical solutions.

# add your code here

Runge-Kutta \(3^{rd}\) order method

Exercise 10.5: Implement the Runge-Kutta \( 3^{\mathrm{rd}} \) order method.

def runge_kutta_3(f, x_0, x_n, y_0, n):

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(runge_kutta_3(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.6: Solve the decay equation for Ra-226 using the Runge-Kutta \( 3^{\mathrm{rd}} \) order method and compare the numerical and analytical solutions.

# add your code here

Runge-Kutta \(4^{th}\) order method

Exercise 10.7: Implement the Runge-Kutta \( 4^{\mathrm{th}} \) order method.

def runge_kutta_4(f, x_0, x_n, y_0, n):

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(runge_kutta_4(f, 0, 10, 2, 10)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.8: Solve the decay equation for Ra-226 using the Runge-Kutta \( 4^{\mathrm{th}} \) order method and compare the numerical and analytical solutions.

# add your code here

Backward (implicit) Euler method

Exercise 10.9: Implement the backward Euler method.

def backward_euler(f, x_0, x_n, y_0, n): # implicit method

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(backward_euler(f, 0, 10, 2, 1000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.10: Solve the equation

\[ \frac{dx}{dy} = - 100 y + 100, \quad y(0) = 5 \]

using the backward Euler method. Compare the numerical solution with the analytical solution

\[ y(x) = (y_0 - 1) e^{-100x} + 1 \]

and with the solution obtained by the forward Euler method.

# add your code here

Multistep methods

Leap-frog method

Exercise 10.11: Implement the Leap-frog method.

def leap_frog(f, x_0, x_n, y_0, n): 

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(leap_frog(f, 0, 10, 2, 10000)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.12: Solve the decay equation for Ra-226 using the Leap-frog method and compare the numerical and analytical solutions.

# add your code here

Adams-Bashforth method

The Adams-Bashforth method

Exercise 10.13: Implement the Adams-Bashforth method.

def adams_bashforth(f, x_0, x_n, y_0, n): 

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(adams_bashforth(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.14: Solve the decay equation for Ra-226 using the Adams-Bashforth method and compare the numerical and analytical solutions.

# add your code here

Adams-Moulton method

The Adams-Moulton method

Exercise 10.15: Implement the Adams-Moulton method.

def adams_moulton(f, x_0, x_n, y_0, n): 

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(adams_moulton(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.16: Solve the decay equation for Ra-226 using the Adams-Moulton method and compare the numerical and analytical solutions.

# add your code here

Predictor-corrector methods

The predictor-corrector method

Exercise 10.17: Implement the predictor-corrector method.

def predictor_corrector(f, x_0, x_n, y_0, n):

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(predictor_corrector(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.18: Solve the decay equation for Ra-226 using the predictor-corrector method and compare the numerical and analytical solutions.

# add your code here

Bulirsch-Stoer method

The Bulirsch-Stoer algorithm

Exercise 10.19: Implement the Bulirsch-Stoer method.

def bulirsch_stoer(f, x_0, x_n, y_0, n, m=5): 

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_ivp function.

def f(t, y): 
    return -0.5 * y
try:
    np.testing.assert_almost_equal(bulirsch_stoer(f, 0, 10, 2, 100)[1][-1], integrate.solve_ivp(f, [0, 10], [2]).y[-1][-1], \
        decimal=4)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 10.20: Solve the decay equation for Ra-226 using the Bulirsch-Stoer method and compare the numerical and analytical solutions.

# add your code here

previous

Tutorial 09 - Quadrature

next

Tutorial 11 - Boundary value problems