Heat Equation Example — Method of Lines¶
This example solves the 1-D heat (diffusion) equation
\[\frac{\partial u}{\partial t} = \alpha\,\frac{\partial^2 u}{\partial x^2},
\quad \alpha = 0.01, \quad x \in [0, 1]\]
with Dirichlet boundary conditions \(u(0,t)=u(1,t)=0\) and the initial condition \(u(x,0)=\sin(\pi x)\).
The spatial dimension is discretized with \(N=49\) interior finite-difference
nodes (method of lines), and the resulting stiff ODE system is integrated in time
with differentialequations.ivpag() using the 'BDF' method.
The analytical solution is:
\[u(x, t) = \sin(\pi x)\,e^{-\alpha\pi^2 t}\]
Example Code¶
"""IMSL heat equation example: 1-D heat equation via method of lines.
Solves the 1-D heat (diffusion) equation
∂u/∂t = α ∂²u/∂x², α = 0.01
on x ∈ [0, 1] with:
- IC: u(x, 0) = sin(πx)
- BC: u(0, t) = u(1, t) = 0 (Dirichlet)
Space is discretized with N interior nodes (method of lines) and time
integration is performed with ivpag() using BDF for the resulting stiff ODE.
Analytical solution: u(x, t) = sin(πx) * exp(-α π² t)
Outputs:
- Maximum absolute error vs analytical at t_final printed to stdout
- SVG plot saved to test_output/example_imsl_heat_equation.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from differentialequations import ivpag
_ALPHA = 0.01
_N = 49 # interior grid points
def run_demo_imsl_heat_equation() -> Dict[str, object]:
"""Run IMSL heat equation example using method of lines with ivpag.
Discretizes ∂u/∂t = α ∂²u/∂x² in space with N interior nodes and uses
ivpag (BDF method) for time integration. Compares numerics against the
analytical solution u(x,t) = sin(πx) * exp(-α π² t).
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``x`` (ndarray), ``t_plot``
(list), ``u_numerical`` (list of ndarray), ``max_error`` (float),
and ``plot_path`` (str).
"""
dx = 1.0 / (_N + 1)
x = np.linspace(dx, 1.0 - dx, _N)
alpha_dx2 = _ALPHA / dx ** 2
u_init = np.sin(np.pi * x)
def heat_mol(t: float, u: np.ndarray) -> np.ndarray:
"""Method-of-lines spatial discretization of the heat equation.
Args:
t (float): Time (unused, autonomous in space terms).
u (np.ndarray): Interior node temperatures of shape ``(N,)``.
Returns:
np.ndarray: Time derivatives du/dt of shape ``(N,)``.
"""
dudt = np.empty_like(u)
dudt[0] = alpha_dx2 * (-2.0 * u[0] + u[1])
dudt[-1] = alpha_dx2 * (u[-2] - 2.0 * u[-1])
dudt[1:-1] = alpha_dx2 * (u[:-2] - 2.0 * u[1:-1] + u[2:])
return dudt
t_plot_targets = [0.1, 0.5, 1.0]
t_final = t_plot_targets[-1]
result = ivpag(heat_mol, (0.0, t_final), u_init, method="BDF", rtol=1e-7, atol=1e-10)
def analytical(x_arr: np.ndarray, t: float) -> np.ndarray:
"""Analytical heat equation solution for the sine initial condition.
Args:
x_arr (np.ndarray): Spatial positions.
t (float): Time.
Returns:
np.ndarray: u(x, t) = sin(πx) * exp(-α π² t).
"""
return np.sin(np.pi * x_arr) * np.exp(-_ALPHA * np.pi ** 2 * t)
u_at_times = []
for t_target in t_plot_targets:
idx = np.argmin(np.abs(result.t - t_target))
u_at_times.append(result.y[:, idx])
u_final_num = result.y[:, -1]
u_final_ana = analytical(x, t_final)
max_error = float(np.max(np.abs(u_final_num - u_final_ana)))
print("\nIMSL Heat Equation Example: Method of Lines with ivpag (BDF)")
print("=" * 60)
print(f" α (diffusivity) : {_ALPHA}")
print(f" Spatial nodes N : {_N} (Δx = {dx:.4f})")
print(f" t_final : {t_final}")
print(f" IC : u(x,0) = sin(πx)")
print(f" BC : u(0,t) = u(1,t) = 0")
print()
print(f" Steps taken : {result.n_steps}")
print(f" Function evals : {result.n_eval}")
print()
print(f" Max |error| vs analytical at t={t_final}: {max_error:.4e}")
print("=" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_heat_equation.svg"
x_full = np.linspace(0, 1, 200)
colors = ["#0369a1", "#0891b2", "#16a34a", "#e11d48"]
t_labels = [0.0] + t_plot_targets
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left panel: solution at multiple times
axes[0].plot(x_full, np.sin(np.pi * x_full),
color=colors[0], linewidth=2.0, label="t = 0 (IC)")
for (t_val, u_num, col) in zip(t_plot_targets, u_at_times, colors[1:]):
u_ana = analytical(x_full, t_val)
axes[0].plot(x_full, u_ana, color=col, linewidth=1.5, linestyle="--", alpha=0.7)
axes[0].plot(x, u_num, color=col, linewidth=2.0, label=f"t = {t_val}")
axes[0].set_xlabel("x")
axes[0].set_ylabel("u(x, t)")
axes[0].set_title("1-D Heat Equation: u(x,t) at selected times\n(dashed = analytical)")
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)
# Right panel: absolute error at t_final
error = np.abs(u_final_num - u_final_ana)
axes[1].semilogy(x, error, color="#e11d48", linewidth=2.0,
label=f"t = {t_final}")
axes[1].set_xlabel("x")
axes[1].set_ylabel(r"$|u_\mathrm{num} - u_\mathrm{ana}|$")
axes[1].set_title(f"Absolute Error vs Analytical (t = {t_final})\nmax error = {max_error:.2e}")
axes[1].legend(fontsize=9)
axes[1].grid(True, which="both", alpha=0.3)
fig.suptitle(
r"IMSL ivpag (BDF): 1-D Heat Equation $\partial_t u = \alpha\,\partial_{xx}u$",
fontsize=13, fontweight="bold"
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"x": x,
"t_plot": t_plot_targets,
"u_numerical": u_at_times,
"max_error": max_error,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_heat_equation()
Plot Output¶
Left: numerical solution (solid) and analytical (dashed) at \(t = 0, 0.1, 0.5, 1.0\). Right: absolute pointwise error versus the analytical solution at \(t = 1.0\).¶
Console Output¶
IMSL Heat Equation Example: Method of Lines with ivpag (BDF)
============================================================
╬▒ (diffusivity) : 0.01
Spatial nodes N : 49 (Δx = 0.0200)
t_final : 1.0
IC : u(x,0) = sin(πx)
BC : u(0,t) = u(1,t) = 0
Steps taken : 19
Function evals : 38
Max |error| vs analytical at t=1.0: 2.9215e-05
============================================================