MATLAB Code#
The MATLAB code used in this book is given here for reference.
Initial value problem solver#
The solveIVP() 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. The functions return arrays containing the values of \(\vec{t}\) and \(\vec{y}\).
function [t, y] = solveIVP(f, tspan, y0, h, solver)
% Define t and y arrays
t = (tspan(1) : h : tspan(2));
y = zeros(length(t), length(y0));
y(1,:) = y0;
% Loop through the steps and calculate single step solver solution
for n = 1 : length(t) - 1
y(n+1,:) = solver(f, t(n), y(n,:), h);
end
end
Explicit Methods#
The Euler method#
function ynew = euler(f, t, y, h)
ynew = y + h * f(t, y);
end
Heun’s method#
function ynew = heun(f, t, y, h)
k1 = f(t, y);
k2 = f(t + h, y + h * k1);
ynew = y + h / 2 * (k1 + k2);
end
The RK4 method#
function ynew = rk4(f, t, y, h)
k1 = f(t, y);
k2 = f(t + 0.5 * h, y + 0.5 * h * k1);
k3 = f(t + 0.5 * h, y + 0.5 * h * k2);
k4 = f(t + h, y + h * k3);
ynew = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
end
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}\).
function [t, y] = solveIVP_SSC(f, tspan, y0, h, tol, solver)
% Define t and y arrays
t = zeros(10000);
y = zeros(10000, length(y0));
t(1) = tspan(1);
y(1,:) = y0;
% Loop through the steps
n = 1;
while t(n) < tspan(2)
% Calculate order p and p+1 solutions
[yp1, yp, order] = solver(f, t(n), y(n,:), h);
% Calculate error estimate
delta = max(abs(yp1 - yp));
% Determine whether the step was successful or not
if delta < tol
y(n+1,:) = yp1;
t(n+1) = t(n) + h;
n = n + 1;
end
% Calculate new value of h (making sure not to exceed tmax)
h = h * max(0.5, min(2, 0.9 * (tol / delta) ^ (1 / (order + 1))));
h = min(h, tspan(2) - t(n));
end
% Remove unused entries from t and y
t(n+1:end) = [];
y(n+1:end,:) = [];
end
Fehlberg’s RKF4(5) method#
function [y5, y4, p, s] = rkf45(f, t, y, h)
k1 = f(t, y);
k2 = f(t + 1/4 * h, y + 1/4 * h * k1);
k3 = f(t + 3/8 * h, y + h * (3/32 * k1 + 9/32 * k2));
k4 = f(t + 12/13 * h, y + h * (1932/2197 * k1 - 7200/2197 * k2 + 7296/2197 * k3));
k5 = f(t + h, y + h * (439/216 * k1 - 8 * k2 + 3680/513 * k3 - 845/4104 * k4));
k6 = f(t + 1/2 * h, y + h * (-8/27 * k1 + 2 * k2 - 3544/2565 * k3 + 1859/4104 * k4 - 11/40 * k5));
y5 = y + h * (16/135 * k1 + 6656/12825 * k3 + 28561/56430 * k4 - 9/50 * k5 + 2/55 * k6);
y4 = y + h * (25/216 * k1 + 1408/2565 * k3 + 2197/4104 * k4 - 1/5 * k5);
p = 4;
s = 6;
end
Defining and solving an IVP#
The following code defines the following ODE and uses the solveIVP() and euler() functions with a step length of \(h=0.2\) to compute the solution
% Define ODE function
f = @(t, y) t * y;
% Define IVP parameters
tspan = [0, 1]; % boundaries of the t domain
y0 = 1; % initial value of the solution
h = 0.2; % step length
% Calculate the solution to the IVP
[t, y] = solveIVP(f, tspan, y0, h, @euler);
Plotting the solution#
The following code uses matplotlib functions to plot the solution.
plot(t, y, 'b-o', LineWidth=1, MarkerFaceColor='b')
axis padded
xlabel('$t$', FontSize=14, Interpreter='latex')
ylabel('$y$', FontSize=14, Interpreter='latex')
legend('Exact', 'RK4', Location='northwest', FontSize=12)
Implicit Runge-Kutta Methods#
Third-order Radau IA IRK method#
function ynew = radauIA(f, t, y, h, tol)
neq = length(y);
Y1 = ones(neq);
Y2 = ones(neq);
Y1old = ones(neq);
Y2old = ones(neq);
for k = 1 : 10
Y1 = y + h * (1/4 * f(t, Y1) - 1/4 * f(t + 2/3 * h, Y2));
Y2 = y + h * (1/4 * f(t, Y1) + 5/12 * f(t + 2/3 * h, Y2));
if max(max(abs(Y1 - Y1old)), max(abs(Y2 - Y2old))) < 1e-6
break
end
Y1old = Y1;
Y2old = Y2;
end
ynew = y + h / 4 * (f(t, Y1) + 3 * f(t + 2/3 * h, Y2));
end
Using MATLAB to solve order conditions#
The following code uses SymPy commands to solve the following order conditions where \(c_2 = 1\)
% Declare symbolic variables
syms a21 b1 b2 c2
% Define known coefficients
c2 = 1;
% Define order conditions
eq1 = b1 + b2 == 1;
eq2 = c2 * b2 == 1/2;
% Define row-sum conditions
eq3 = a21 == c2;
% Solve order conditions
solve(eq1, eq2, eq3)
Stability#
Plotting stability regions#
The following code plots the region of absolute stability for the Euler method.
% Generate z values
xmin = -3;
xmax = 1;
ymin = -1.5;
ymax = 1.5;
[X, Y] = meshgrid(linspace(xmin, xmax, 200), linspace(ymin, ymax, 200));
Z = X + Y * 1i;
% Define stability function
R = 1 + Z;
% Plot stability region
contourf(X, Y, abs(R), [0, 1], LineWidth=2)
xline(0, LineWidth=2)
yline(0, LineWidth=2)
colormap([153, 204, 255 ; 255, 255, 255] / 255)
axis equal
axis([xmin, xmax, ymin, ymax])
xlabel("$\mathrm{Re}(z)$", FontSize=12, Interpreter="latex")
ylabel("$\mathrm{Im}(z)$", FontSize=12, Interpreter="latex")
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 = [0, 0, 0, 0 ;
1/2, 0, 0, 0 ;
0, 1/2, 0, 0 ;
0, 0, 1, 0];
b = [1/6 ; 1/3 ; 1/3 ; 1/6];
e = ones(4, 1);
% Calculate coefficients
for k = 1 : 5
sym(b' * A ^ (k - 1) * e)
end
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 IRK method
A = [5/12, -1/12 ; 3/4, 1/4];
ebT = [3/4, 0 ; 0, 1/4];
I = eye(size(A, 1));
% Define numerator and denominator functions
P = @(z) det(I - z * A + z * ebT);
Q = @(z) det(I - z * A);
% Calculate R(z)
syms z y
Rz = P(z) / Q(z)
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_of_Q = solve(Q(z) == 0)
% Check E(y) >= 0
E = Q(1i * y) * Q(-1i * y) - P(1i * y) * P(-1i * y);
E = simplify(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\).
function [L, U] = lu(A)
n = size(A, 1);
L = eye(n);
U = zeros(n);
for j = 1 : n
for i = 1 : n
sum_ = 0;
if i <= j
for k = 1 : i - 1
sum_ = sum_ + L(i,k) * U(k,j);
end
U(i,j) = A(i,j) - sum_;
else
for k = 1 : j - 1
sum_ = sum_ + L(i,k) * U(k,j);
end
L(i,j) = (A(i,j) - sum_) / U(j,j);
end
end
end
end
Forward and back substitution#
The following code defines the functions forward_substitution() and back_substitution() which perform forward and back substitution.
function x = forward_substitution(L, b)
n = size(L, 1);
x = zeros(n, 1);
for i = 1 : n
sum_ = 0;
for j = 1 : i - 1
sum_ = sum_ + L(i,j) * x(j);
end
x(i) = (b(i) - sum_) / L(i,i);
end
end
function x = back_substitution(U, b)
n = size(U, 1);
x = zeros(n, 1);
x(n) = b(n) / U(n,n);
for i = n - 1 : -1 : 1
sum_ = 0
for j = i + 1 : n
sum_ = sum_ + U(i,j) * x(j);
end
x(i) = (b(i) - sum_) / U(i,i);
end
end
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.
function P = partial_pivot(A)
n = size(A, 1);
P = eye(n);
for j = 1 : n
maxpivot = abs(A(j,j));
maxpivotrrow = j;
for i = j + 1 : n
if abs(A(i,j)) > maxpivot
maxpivot = abs(A(i,j));
maxpivotrow = i;
end
end
A([j,maxpivotrow],:) = A([maxpivotrow,j],:);
P([j,maxpivotrow],:) = P([maxpivotrow,j],:);
end
end
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}\).
function L = cholesky(A)
n = size(A, 1);
L = zeros(n);
for j = 1 : n
for i = j : n
sum_ = 0;
for k = 1 : j - 1
sum_ = sum_ + L(i,k) * L(j,k);
end
if i == j
L(i,j) = sqrt(A(i,j) - sum_);
else
L(i,j) = (A(i,j) - sum_) / L(j,j);
end
end
end
end
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\).
function [Q, R] = qr_gramschmidt(A)
n = size(A, 2);
Q = zeros(size(A));
R = zeros(n);
for j = 1 : n
sum_ = 0;
for i = 1 : j - 1
R(i,j) = dot(Q(:,i), A(:,j));
sum_ = sum_ + R(i,j) * Q(:,i);
end
u = A(:,j) - sum_;
R(j,j) = norm(u);
Q(:,j) = u / R(j,j);
end
end
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\).
function [Q, R] = qr_householder(A)
[m, n] = size(A);
I = eye(m);
Q = I;
R = A;
for j = 1 : n
u = R(:,j);
u(1:j-1) = 0;
u = u + sign(R(j,j)) * norm(u) * I(:,j);
v = u / norm(u);
H = I - 2 * v * v';
Q = Q * H;
R = H * R;
end
end
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.
function lambda = eigvals(A, tol)
for k = 1 : 20
[Q, R] = qr_householder(A);
Aprev = A;
A = R * Q;
if max(abs(diag(A - Aprev))) < tol
break
end
end
lambda = diag(A);
end
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#
function x = jacobi(A, b, tol)
n = length(b);
x = zeros(n, 1);
xnew = zeros(n, 1);
for k = 1 : 100
for i = 1 : n
sum_ = 0;
for j = 1 : n
if j ~= i
sum_ = sum_ + A(i,j) * x(j);
end
end
xnew(i) = (b(i) - sum_) / A(i,i);
end
x = xnew;
r = b - A * x;
if max(abs(r)) < tol
break
end
end
end
The Gauss-Seidel method#
function x = gauss_seidel(A, b, tol)
n = length(b);
x = zeros(n, 1);
for k = 1 : 100
for i = 1 : n
sum_ = 0;
for j = 1 : n
if j ~= i
sum_ = sum_ + A(i,j) * x(j);
end
end
x(i) = (b(i) - sum_) / A(i,i);
end
r = b - A * x;
if max(abs(r)) < tol
break
end
end
end
The SOR method#
function x = sor(A, b, omega, tol)
n = length(b);
x = zeros(n, 1);
for k = 1 : 100
for i = 1 : n
sum_ = 0;
for j = 1 : n
if j ~= i
sum_ = sum_ + A(i,j) * x(j);
end
end
x(i) = (1 - omega) * x(i) + omega * (b(i) - sum_) / A(i,i);
end
r = b - A * x;
if max(abs(r)) < tol
break
end
end
end