Python Code#
The Python code used in this book is given here for reference.
Importing libraries#
The following code is used to import the libraries that we will use here.
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
NumPy (pronounced ‘num-pie’) for performing numerical calculations
SymPy (pronounced ‘sim-pie’) for performing symbolic calculations
matplotlib for plotting solutions
Explicit Methods#
The Euler method#
def euler(f, tspan, y0, h):
nsteps = int((tspan[1] - tspan[0]) / h)
N = len(y0)
t = np.zeros(nsteps + 1)
y = np.zeros((N, nsteps + 1))
t[0] = tspan[0]
y[:,0] = y0
for n in range(nsteps):
y[:,n+1] = y[:,n] + h * f(t[n], y[:,n])
t[n+1] = t[n] + h
return t, y.T
The RK4 method#
def rk4(f, tspan, y0, h):
nsteps = int((tspan[1] - tspan[0]) / h)
N = len(y0)
t = np.zeros(nsteps + 1)
y = np.zeros((N, nsteps + 1))
t[0] = tspan[0]
y[:,0] = y0
for n in range(nsteps):
k1 = f(t[n], y[:,n])
k2 = f(t[n] + 0.5 * h, y[:,n] + 0.5 * h * k1)
k3 = f(t[n] + 0.5 * h, y[:,n] + 0.5 * h * k2)
k4 = f(t[n] + h, y[:,n] + h * k3)
y[:,n+1] = y[:,n] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
t[n+1] = t[n] + h
return t, y.T
Adaptive step size control#
The solveIVP_SSC() function is used to solve an initial value problem of the form \(\vec{y}' = f(t, \vec{y})\), \(t \in [t_{\min}, t_{\max}]\) and \(\vec{y}_0 = \vec{\alpha}\) using a single step method with adaptive step size control. The functions return arrays containing the values of \(\vec{t}\) and \(\vec{y}\).
Bogacki-Shampine 3(2) Runge-Kutta method#
def rk23(f, tspan, y0, atol=1e-6, rtol=1e-3):
N = len(y0)
t = np.zeros(100000)
y = np.zeros((N, 100000))
t[0] = tspan[0]
y[:,0] = y0
h = 0.8 * rtol ** (1 / 3)
n = 0
while t[n] < tspan[-1]:
k1 = f(t[0], y[:,0])
k2 = f(t[n] + 1/2 * h, y[:,n] + 1/2 * h * k1)
k3 = f(t[n] + 3/4 * h, y[:,n] + 3/4 * h * k2)
k4 = f(t[n] + h, y[:,n] + h * (2/9 * k1 + 1/3 * k2 + 4/9 * k3))
y3 = y[:,n] + h * (2/9 * k1 + 1/3 * k2 + 4/9 * k3)
y2 = y[:,n] + h * (7/24 * k1 + 1/4 * k2 + 1/3 * k3 + 1/8 * k4)
delta = np.linalg.norm(y3 - y2)
tol = atol + rtol * np.linalg.norm(y3)
if delta < tol:
y[:,n+1] = y3
t[n+1] = t[n] + h
k1 = k4
n += 1
r = max(0.1, 0.8 * (tol / delta) ** (1/3))
h = min(r * h, tspan[-1] - t[n])
return t[:n+1], y[:,:n+1].T
Plotting the solution#
The following code uses matplotlib functions to plot the solution.
fig, ax = plt.subplots()
plt.plot(t, y, "bo-", label="Euler")
plt.xlabel("$t$")
plt.ylabel("$y$")
plt.legend()
plt.show()
Implicit Runge-Kutta Methods#
Third-order Radau IA IRK method#
def jac(f, t, y):
J = np.zeros((len(y), len(y)))
epsilon = 1e-6
for i in range(len(y)):
y_plus_epsilon = y.astype(float)
y_minus_epsilon = y.astype(float)
y_plus_epsilon[i] += epsilon
y_minus_epsilon[i] -= epsilon
J[:,i] = (f(t, y_plus_epsilon) - f(t,y_minus_epsilon)) / (2 * epsilon)
return J
def radauIA(f, tpsan, y0, h):
N = len(y0)
nsteps = int((tspan[1] - tspan[0]) / h)
t = np.zeros(nsteps + 1)
y = np.zeros((N, nsteps + 1))
t[0] = tspan[0]
y[:,0] = y0
A = np.array([[1/4, -1/4],
[1/4, 5/12]])
b = np.array([1/4, 3/4])
c = np.array([0, 2/3])
s = 2
for n in range(nsteps):
z = np.zeros(N * s)
F = np.zeros(N * s)
J = jac(f, t[n], y[:,n])
for k in range(10):
F[:N] = f(t[n] + c[0] * h, y[:,n] + z[:N])
F[N:] = f(t[n] + c[1] * h, y[:,n] + z[N:])
g = z - h * np.dot(np.kron(A, np.eye(N)), F)
delta_z = np.linalg.solve(np.eye(N * s) - h * np.kron(A, J), -g)
z += delta_z
if np.linalg.norm(delta_z) < 1e-6:
break
y[:,n+1] = y[:,n] + h * np.dot(np.kron(b.T, np.eye(N)), F)
t[n+1] = t[n] + h
return t, y.T
Using SymPy to solve order conditions#
The following code uses SymPy commands to solve the following order conditions where \(c_2 = 1\)
import sympy as sp
# Declare symbolic variables
a21, b1, b2, c2 = sp.symbols("a21, b1, b2, c2")
# Define known coefficients
c2 = 1
# Define order conditions
eq1 = b1 + b2 - 1
eq2 = b2 * c2 - sp.Rational(1,2)
# Define row sum conditions
eq3 = a21 - c2
# Solve order conditions
sp.solve((eq1, eq2, eq3), manual=True)
Stability#
Plotting stability regions#
The following code plots the region of absolute stability for the Euler method.
# Generate z values
xmin, xmax, ymin, ymax = -3, 1, -1.5, 1.5
X, Y = np.meshgrid(np.linspace(xmin, xmax, 200),np.linspace(ymin, ymax, 200))
Z = X + Y * 1j
# Define stability function
R = 1 + Z
# Plot stability region
fig = plt.figure()
contour = plt.contourf(X, Y, abs(R), levels=[0, 1], colors="#99ccff")
plt.contour(X, Y, abs(R), colors= "k", levels=[0, 1])
plt.axhline(0, color="k", linewidth=1)
plt.axvline(0, color="k", linewidth=1)
plt.axis("equal")
plt.axis([xmin, xmax, ymin, ymax])
plt.xlabel("$\mathrm{Re}(z)$", fontsize=12)
plt.ylabel("$\mathrm{Im}(z)$", fontsize=12)
plt.show()
Stability function for explicit methods#
The following code calculates the stability function for an explicit Runge-Kutta method defined by the following Butcher tableau
# Define ERK method
A = sp.Matrix([[0, 0, 0, 0],
[sp.Rational(1,2), 0, 0, 0],
[0, sp.Rational(1,2), 0, 0],
[0, 0, 1, 0]])
b = sp.Matrix([sp.Rational(1,6), sp.Rational(1,3), sp.Rational(1,3), sp.Rational(1,6)])
e = sp.ones(4,1)
# Determine coefficients of the stability function
for k in range(4):
display(b.T * A ** k * e)
Stability function for implicit methods#
The following code calculates the stability function for an explicit Runge-Kutta method defined by the following Butcher tableau.
# Define numerator and denominator functions
def P(z):
return (I - z * A + z * ebT).det()
def Q(z):
return (I - z * A).det()
# Define IRK method
A = sp.Matrix([[sp.Rational(5,12), -sp.Rational(1,12)],
[sp.Rational(3,4), sp.Rational(1,4)]])
ebT = sp.Matrix([[sp.Rational(3,4), 0], [0, sp.Rational(1,4)]])
I = sp.eye(2)
# Calculate R(z)
z, y = sp.symbols('z, y')
R = P(z) / Q(z)
print(f"R(z) = ")
display(R)
Checking if an implicit method is A-stable#
The following code outputs the following conditions for A-stability
The roots of \(Q(z)\) have positive real parts
\(E(y) = Q(iy)Q(-iy) - P(iy)P(-iy) \geq 0\)
where the stability function for the method is \(R(z) = \dfrac{P(z)}{Q(z)}\).
# Check roots of Q have positive real parts
roots = sp.solve(Q(z) - 0)
print(f"Roots of Q(z)")
display(roots)
# Check E(y) >= 0
E = Q(1j * y) * Q(-1j * y) - P(1j * y) * P(-1j * y)
print(f"E(y) = ")
display(sp.simplify(sp.nsimplify(E)))
Matrix Decomposition Methods#
LU decomposition#
The following code defines the function lu() which calculates the LU decomposition of a square matrix \(A\) and returns the lower and upper triangular matrices \(L\) and \(U\) such that \(A = LU\).
def lu(A):
n = A.shape[0]
L, U = np.eye(n), np.zeros((n, n))
for j in range(n):
for i in range(n):
sum_ = 0
if i <= j:
for k in range(i):
sum_ += L[i,k] * U[k,j]
U[i,j] = A[i,j] - sum_
else:
for k in range(j):
sum_ += L[i,k] * U[k,j]
L[i,j] = (A[i,j] - sum_) / U[j,j]
return L, U
Forward and back substitution#
The following code defines the functions forward_substitution() and back_substitution() which perform forward and back substitution.
def forward_substitution(L, b):
n = L.shape[0]
x = np.zeros(n)
for i in range(n):
sum_ = 0
for j in range(i):
sum_ += L[i,j] * x[j]
x[i] = (b[i] - sum_) / L[i,i]
return x
def back_substitution(U, b):
n = U.shape[0]
x = np.zeros(n)
x[-1] = b[-1] / U[-1,-1]
for i in range(n - 2, -1, -1):
sum_ = 0
for j in range(i + 1, n):
sum_ += U[i,j] * x[j]
x[i] = (b[i] - sum_) / U[i,i]
return x
Partial pivoting#
The following code defines the function partial_pivot() that performs partial pivoting on a matrix and outputs the matrix and the permutation matrix.
def partial_pivot(A):
n = A.shape[0]
P = np.eye(n)
for j in range(n):
maxpivot, maxpivotrow = abs(A[j,j]), j
for i in range(j + 1, n):
if abs(A[i,j]) > maxpivot:
maxpivot, maxpivotrow = abs(A[i,j]), i
A[[j,maxpivotrow]] = A[[maxpivotrow,j]]
P[[j,maxpivotrow]] = P[[maxpivotrow,j]]
return P
Cholesky decomposition#
The following code defines the function cholesky() which performs Cholesky decomposition on a matrix \(A\) and outputs the lower triangular matrix \(L\) such that \(A = LL^\mathrm{T}\).
def cholesky(A):
n = A.shape[0]
L = np.zeros((n, n))
for j in range(n):
for i in range(j, n):
for k in range(j):
L[i,j] += L[i,k] * L[j,k]
if i == j:
L[i,j] = np.sqrt(A[i,j] - L[i,j])
else:
L[i,j] = (A[i,j] - L[i,j]) / L[j,j]
return L
QR decomposition using the Gram-Schmidt process#
The following code defines the function qr_gramschmidt() which performs QR decomposition using the Gram-Schmidt process on a matrix \(A\) and outputs the orthogonal matrix \(Q\) and upper triangular matrix \(R\) such that \(A = QR\).
def qr_gramschmidt(A):
n = A.shape[1]
Q, R = np.zeros(A.shape), np.zeros((n,n))
for j in range(n):
sum_ = 0
for i in range(j):
R[i,j] = np.dot(Q[:,i], A[:,j])
sum_ += R[i,j] * Q[:,i]
u = A[:,j] - sum_
R[j,j] = np.linalg.norm(u)
Q[:,j] = u / R[j,j]
return Q, R
QR decomposition using the Householder transformations#
The following code defines the function qr_householder() which performs QR decomposition using Household transformations on a matrix \(A\) and outputs the orthogonal matrix \(Q\) and upper triangular matrix \(R\) such that \(A = QR\).
def qr_householder(A):
m, n = A.shape
I = np.eye(m)
Q, R = np.eye(m), np.copy(A)
for j in range(n):
u = R[:,[j]]
u[:j] = 0
u = u + np.sign(R[j,j]) * np.linalg.norm(u) * I[:,[j]]
v = u / np.linalg.norm(u)
H = I - 2 * np.dot(v, v.T)
R = np.dot(H, R)
Q = np.dot(Q, H)
return Q, R
Calculating eigenvalues of a matrix using the QR algorithm.#
The following code defines the function eigvals() which calculates the eigenvalues of a square matrix \(A\) using the QR algorithm.
def eigvals(A, tol=1e-6):
for k in range(20):
Q, R = qr_householder(A)
A, Aprev = np.matmul(R, Q), A
if max(abs(np.diagonal(A - Aprev))) < tol:
break
return np.diagonal(A)
Indirect methods#
The following methods calculate the solutions to the system of linear equations \(A \vec{x} = \vec{b}\) ceasing iterations when the largest value of the residual is less than tol.
The Jacobi method#
def jacobi(A, b, tol=1e-6):
n = len(b)
x, xnew = np.zeros(n), np.zeros(n)
for k in range(100):
for i in range(n):
sum_ = 0
for j in range(n):
if i != j:
sum_ += A[i,j] * x[j]
xnew[i] = (b[i] - sum_) / A[i,i]
x = np.copy(xnew)
r = b - np.dot(A, x)
if max(abs(r)) < tol:
break
return x
The Gauss-Seidel method#
def gauss_seidel(A, b, tol=1e-6):
n = len(b)
x = np.zeros(n)
for k in range(100):
for i in range(n):
sum_ = 0
for j in range(n):
if i != j:
sum_ += A[i,j] * x[j]
x[i] = (b[i] - sum_) / A[i,i]
r = b - np.dot(A, x)
if max(abs(r)) < tol:
break
return x
The SOR method#
def sor(A, b, omega, tol=1e-6):
n = len(b)
x = np.zeros(n)
for k in range(100):
for i in range(n):
sum_ = 0
for j in range(n):
if i != j:
sum_ += A[i,j] * x[j]
x[i] = (1 - omega) * x[i] + omega * (b[i] - sum_) / A[i,i]
r = b - np.dot(A, x)
if max(abs(r)) < tol:
break
return x