Workshop 2: Convergence Analysis of IVP problems#
In this workshop, we will:
verify the fifth-order convergence numerically for the RK45 method.
Convergence analysis for the forward Euler scheme.
introduce adaptive timestepping through RK45.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
Problem definition#
We solve the following ODE -
\[m\ddot{u} + c\dot{u} + ku = 0\]
With the state vector being -
\[\begin{split}
\boldsymbol{q} =
\begin{bmatrix}
u \\
\dot{u}
\end{bmatrix}\end{split}\]
We start by defining the numerical parameters for this mass-spring-damper system.
# System parameters
m = 1.0 # mass [kg]
c = 0.2 # damping [Ns/m]
k = 4.0 # stiffness [N/m]
# Time interval
t_0 = 0.0
t_f = 10.0
# Initial conditions
init_disp = 1.0
init_velo = 0.0
q0 = [init_disp, init_velo]
# Define the system of equations
def q_dot(t, q):
u = q[0]
v = q[1]
du_dt = v
dv_dt = -(k/m) * u - (c/m) * v
return np.array([du_dt, dv_dt])
Exact solution#
We define -
\[
\omega_n = \sqrt{\frac{k}{m}}, \qquad
\zeta = \frac{c}{2\sqrt{km}}.
\]
For the underdamped case \((\zeta < 1)\), the damped frequency is
\[
\omega_d = \omega_n\sqrt{1-\zeta^2}.
\]
The exact displacement is
\[
u(t)=e^{-\zeta\omega_n t}
\left[
u_0\cos(\omega_d t)
+
\frac{v_0+\zeta\omega_n u_0}{\omega_d}\sin(\omega_d t)
\right].
\]
This exact solution allows us to measure the numerical error directly.
We define the exact solution#
omega_n = np.sqrt(k / m)
zeta = c / (2 * np.sqrt(k * m))
if zeta >= 1.0:
raise ValueError("This notebook uses the underdamped exact solution, so choose parameters with zeta < 1.")
omega_d = omega_n * np.sqrt(1 - zeta**2)
def u_exact(t):
A = q0[0]
B = (q0[1] + zeta * omega_n * q0[0]) / omega_d
t = np.asarray(t)
return np.exp(-zeta * omega_n * t) * (
A * np.cos(omega_d * t) + B * np.sin(omega_d * t)
)
def v_exact(t):
A = q0[0]
B = (q0[1] + zeta * omega_n * q0[0]) / omega_d
t = np.asarray(t)
exp_term = np.exp(-zeta * omega_n * t)
return exp_term * (
-zeta * omega_n * (A * np.cos(omega_d * t) + B * np.sin(omega_d * t))
-A * omega_d * np.sin(omega_d * t)
+B * omega_d * np.cos(omega_d * t)
)
# Exact solution state function
def q_exact(t):
return np.vstack((u_exact(t), v_exact(t)))
Convergence analysis#
For a fifth-order Runge–Kutta method, the global error is expected to behave like
\[
\epsilon(\Delta t) \sim C (\Delta t)^5.
\]
To verify this numerically, we solve the same problem on a sequence of timestep grids and compute the error with respect to the exact solution.
We estimate the observed order of convergence by measuring the slope of the best-fit line in a log-log plot of \(\epsilon\) versus \(\Delta t\).
dt_values = np.array([0.1, 0.05, 0.025, 0.0125, 0.00625])
errors_final_time = []
errors_max = []
def solve_at_dt(dt):
t_eval = np.arange(t_0, t_f + 0.5 * dt, dt)
sol = solve_ivp(
fun=q_dot,
t_span=(t_0, t_f),
y0=q0,
t_eval=t_eval,
method='RK45',
max_step=dt # Ensure that the solver does not take steps larger than dt
)
return sol, t_eval
# Loop over different time step sizes and compute errors using the function solve_at_dt
# Write your code here to compute the errors and store them in errors_final_time and errors_max lists
We now plot the errors we have obtained
# log-log fit
if len(errors_final_time) == 0 or len(errors_max) == 0:
raise ValueError(
"errors_final_time and errors_max are empty. Run the previous exercise cell "
"to compute the errors before fitting and plotting."
)
if len(errors_final_time) != len(dt_values) or len(errors_max) != len(dt_values):
raise ValueError(
"errors_final_time and errors_max must have the same length as dt_values."
)
if np.any(np.asarray(errors_final_time) <= 0) or np.any(np.asarray(errors_max) <= 0):
raise ValueError(
"All error values must be strictly positive before taking logarithms."
)
slope_final_timestep, _ = np.polyfit(np.log(dt_values), np.log(errors_final_time), 1)
slope_max, _ = np.polyfit(np.log(dt_values), np.log(errors_max), 1)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[6], line 3
1 # log-log fit
2 if len(errors_final_time) == 0 or len(errors_max) == 0:
----> 3 raise ValueError(
4 "errors_final_time and errors_max are empty. Run the previous exercise cell "
5 "to compute the errors before fitting and plotting."
6 )
ValueError: errors_final_time and errors_max are empty. Run the previous exercise cell to compute the errors before fitting and plotting.
plt.figure(figsize=(6,4))
# data
plt.loglog(dt_values, errors_final_time, "o-", label="Final timestep error")
plt.loglog(dt_values, errors_max, "s-", label="Maximum error")
# fitted lines
plt.loglog(
dt_values,
np.exp(np.polyfit(np.log(dt_values), np.log(errors_final_time), 1)[1]) * dt_values**slope_final_timestep,
"--",
label=f"Best Fit final timestep (p={slope_final_timestep:.2f})"
)
plt.loglog(
dt_values,
np.exp(np.polyfit(np.log(dt_values), np.log(errors_max), 1)[1]) * dt_values**slope_max,
"--",
label=f"Best Fit maximum error (p={slope_max:.2f})"
)
plt.xlabel(r"Timestep $\Delta t$")
plt.ylabel(r"Error $\epsilon$")
plt.title("Observed order of convergence (RK5)")
plt.legend()
plt.show()
We now perform the convergence analysis for the Forward Euler method using the FE_solver implementation we observed in Workshop 1.
# FE solver
def FE_solver(qdot, tspan, q0):
# Initialize the solution vector at the times we want to evaluate.
# Match solve_ivp's state layout: one column per timestep, one row per state variable.
nsteps = len(tspan)
q = np.zeros((2, nsteps))
# Fill initial step
q[:, 0] = q0
# Loop over time entries
for i in range(nsteps - 1):
# Evaluate the derivative using the function provided as argument
dt_step = tspan[i + 1] - tspan[i]
dq = qdot(tspan[i], q[:, i])
# Update the solution vector at step "i"
q[:, i + 1] = q[:, i] + np.multiply(dq, dt_step)
return q
dt_values = np.array([0.1, 0.05, 0.025, 0.0125, 0.00625])
errors_final_time = []
errors_max = []
def solve_at_dt_FE_solver(dt):
t_eval = np.arange(t_0, t_f + 0.5 * dt, dt)
sol = FE_solver(
qdot=q_dot,
tspan=t_eval,
q0=q0
)
return sol, t_eval
# Loop over different time step sizes and compute errors using the function solve_at_dt_FE_solver
## Write your code here to compute the errors and store them in errors_final_time and errors_max lists
# log-log fit
slope_final_timestep, intercept_final = np.polyfit(np.log(dt_values), np.log(errors_final_time), 1)
slope_max, intercept_max = np.polyfit(np.log(dt_values), np.log(errors_max), 1)
plt.figure(figsize=(6,4))
# error data
plt.loglog(dt_values, errors_final_time, "o-", label="Final-time error")
plt.loglog(dt_values, errors_max, "s-", label="Maximum error")
# fitted lines
plt.loglog(
dt_values,
np.exp(intercept_final) * np.array(dt_values)**slope_final_timestep,
"--",
label=f"Best fit final-time error (p={slope_final_timestep:.2f})"
)
plt.loglog(
dt_values,
np.exp(intercept_max) * np.array(dt_values)**slope_max,
"--",
label=f"Best fit maximum error (p={slope_max:.2f})"
)
plt.xlabel(r"Timestep $\Delta t$")
plt.ylabel(r"Error $\epsilon$")
plt.title("Observed order of convergence (Forward Euler)")
plt.legend()
plt.show()
Adaptive timestepping#
In the previous section, we enforced a restricted timestep using max_step to study convergence with respect to \(\Delta t\).
In practice, however, solve_ivp with RK45 does not use a fixed timestep. Instead, it automatically adjusts the timestep based on an estimate of the local error. This is known as adaptive timestepping.
# Function to solve the system using an adaptive time-stepping method (RK45)
def solve_adaptive(rtol=1e-6, atol=1e-9):
sol = solve_ivp(
fun=q_dot,
t_span=(t_0, t_f),
y0=q0,
method="RK45",
rtol=rtol,
atol=atol
)
return sol
We now vary the relative tolerance in solve_adaptive to investigate the number of timesteps selected by the solver for different tolerance values selected.
# Tolerance values to test
rtol_values = [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]
nsteps = []
errors_final_time = []
errors_max = []
for rtol in rtol_values:
atol = rtol * 1e-3
sol_i = solve_adaptive(rtol=rtol, atol=atol)
# Exact solution evaluated at the adaptive time points
q_ex = q_exact(sol_i.t)
# Pointwise error in time
error_time = np.linalg.norm(sol_i.y - q_ex, axis=0)
# Number of accepted timesteps
nsteps.append(len(sol_i.t) - 1)
# Error measures
errors_final_time.append(error_time[-1])
errors_max.append(np.max(error_time))
plt.figure(figsize=(7,4))
plt.loglog(rtol_values, nsteps, "o-")
plt.xlabel("Relative tolerance")
plt.ylabel("Number of timesteps")
plt.title("Adaptive RK45: timesteps versus tolerance")
plt.show()