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
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{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} $$Note: In the case of explicit methods the Runge-Kutta matrix is lower triangular.
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
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
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
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
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
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
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
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
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
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