.ipynb .pdf

repository open issue

Binder Live Code

Contents

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),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),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{array}{r}\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \cdots & a_{1r}\\ c_2& a_{21}& a_{22}& \cdots & a_{2r}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_r& a_{r1}& a_{r2}& \cdots & a_{rr}\\ & b_1& b_2& \cdots & b_r\end{array}\end{array}$$

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),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}=-100y+100,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