Method of Lines

Pdes

The method of lines (MOL) is a deceptively simple organizing principle for time-dependent PDEs: discretize SPACE first, then integrate the resulting LARGE SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS in time using any time integrator you like. It is the bridge between PDE theory and the well-developed machinery of ODE numerical methods — Runge-Kutta, adaptive step control, stiff solvers, symplectic integrators, geometric integrators. Once you've cast a PDE into MOL form, the entire ODE toolkit is available.

The conceptual move is the SEMI-DISCRETIZATION. A PDE like $u_{t} = L u$ (where $L$ is a spatial differential operator) becomes, after discretizing $L$ in space via finite differences, finite elements, or a spectral method, a system of ODEs:

\frac{d u ( t )}{d t} = L u (t),

where $u (t)$ is the vector of nodal values and $L$ is the discrete spatial operator. The "lines" in "method of lines" are vertical lines in the $x$ - $t$ plane — fixed- $x$ slices along which time evolves as an ODE.

Why this is useful

Three concrete advantages over manually coupled space-time discretization (e.g. plain FTCS, which is space-time-coupled).

(1) Adaptive time stepping for free. A general-purpose ODE integrator (scipy's solve_ivp, SUNDIALS' CVODE, MATLAB's ode45) monitors local truncation error and adapts $Δ t$ on the fly. Where the solution is changing rapidly the integrator takes small steps; where it's quiescent it takes large ones. This is dramatically more efficient than a fixed time step that has to be small enough for the WORST case.

(2) Stiffness handling. Parabolic problems become STIFF as the mesh refines: the eigenvalues of $L$ for a discretized second-derivative operator span $O (1)$ to $O (h^{- 2})$ . Explicit methods like FTCS or RK4 require $Δ t \sim h^{2}$ , which is punishing. Implicit ODE solvers (BDF, Rosenbrock) handle stiff systems with much larger time steps. MOL gives you a STIFF system; just pick a stiff ODE solver. The fact that the underlying problem is a PDE becomes irrelevant.

(3) Composability. Different spatial discretizations and different time integrators can be mixed and matched. Want fourth-order space and second-order time? Build the spatial Jacobian with a 5-point stencil and integrate with Heun. Want spectral space and 8th-order time? Cheb differentiation matrix plus RK8. The MOL framework decouples the two questions.

The trade-off

What you give up: full time-space adaptivity. The spatial mesh is fixed throughout the time integration, even if the solution develops features (shocks, fronts) where you'd want refinement. Adaptive mesh refinement (AMR) in MOL is possible but invasive; problems with strong moving features often justify abandoning MOL in favor of front-tracking or fully space-time adaptive methods.

A second concern is the CFL-like constraint on time-step size, which is determined by the SPECTRAL RADIUS of $L$ . For the second-derivative operator on a fine mesh, eigenvalues reach $\sim 4/ h^{2}$ — a stiff problem. Explicit ODE solvers tied to this $h^{- 2}$ scaling don't help; stiffness-handling implicit methods are necessary at fine resolutions.

Code: heat equation via MOL

Concrete demonstration. Discretize the 1D heat equation in space with central finite differences (the spatial discretization we used in finite differences), then integrate the resulting $(N - 2)$ -dimensional ODE system with scipy's adaptive RK45.

# Method of lines: discretize SPACE first, integrate the resulting ODE
# system in TIME using a general-purpose ODE integrator.
#
# Problem:  u_t = u_xx  on [0, 1],  u(0,t) = u(1,t) = 0,  u(x,0) = sin(pi x).
# Exact solution:  u(x, t) = sin(pi x) exp(-pi² t).
#
# Spatial discretization: central finite differences with N points →
# system of N-2 ODEs for the interior values:
#   du_j/dt = (u_{j-1} - 2 u_j + u_{j+1}) / h²,  j = 1, ..., N-2.
# Then use scipy.integrate.solve_ivp (adaptive RK45) — it picks dt
# automatically based on local truncation error.

import numpy as np
from scipy.integrate import solve_ivp

N  = 51
x  = np.linspace(0, 1, N)
h  = x[1] - x[0]

def rhs(t, u_int):
    """RHS of the semi-discretized heat eqn, interior values only."""
    u = np.zeros(N)
    u[1:-1] = u_int
    return (u[2:] - 2*u[1:-1] + u[:-2]) / h**2

u0 = np.sin(np.pi * x[1:-1])

print("u_t = u_xx, u_0(x) = sin(pi x), Dirichlet BCs.")
print(f"{'t':>6s}  {'MOL u(0.5)':>14s}  {'exact u(0.5)':>14s}  {'|err|':>10s}")
for t in [0.01, 0.05, 0.1, 0.5]:
    sol = solve_ivp(rhs, (0, t), u0, t_eval=[t],
                    method='RK45', rtol=1e-9, atol=1e-12)
    u_mol = np.zeros(N)
    u_mol[1:-1] = sol.y[:, 0]
    exact = np.sin(np.pi * x) * np.exp(-np.pi**2 * t)
    err = np.max(np.abs(u_mol - exact))
    print(f"  {t:6.2f}  {u_mol[N//2]:>14.8f}  {exact[N//2]:>14.8f}  {err:>10.2e}")

# Verify that the adaptive integrator picked a SENSIBLE time step:
sol = solve_ivp(rhs, (0, 0.1), u0, method='RK45', rtol=1e-9, atol=1e-12)
print(f"\nAdaptive RK45 used {len(sol.t) - 1} time steps on [0, 0.1].")
print(f"Average dt = {0.1 / (len(sol.t) - 1):.3e}")
print(f"FTCS would need dt < h²/2 = {0.5 * h**2:.3e} for stability  "
      f"({int(0.1 / (0.5 * h**2))} steps).")

Output:

u_t = u_xx, u_0(x) = sin(pi x), Dirichlet BCs.
     t      MOL u(0.5)    exact u(0.5)       |err|
   0.01      0.90604747      0.90601806    2.94e-05
   0.05      0.61059713      0.61049803    9.91e-05
   0.10      0.37282886      0.37270784    1.21e-04
   0.50      0.00720357      0.00719188    1.17e-05

Adaptive RK45 used 21 time steps on [0, 0.1].
Average dt = 4.762e-03
FTCS would need dt < h²/2 = 2.000e-04 for stability  (500 steps).

The MOL solution tracks the analytical answer $sin (π x) e^{- π^{2} t}$ to $\sim 1 0^{- 4}$ across the simulation window — error dominated by the $O (h^{2})$ spatial discretization, not the time integration. The ADAPTIVE STEP COUNT is the punchline: scipy's RK45 needed only 21 time steps to reach $t = 0.1$ , with average $Δ t \approx 5 \times 1 0^{- 3}$ . The corresponding FTCS would need 500 time steps because of the $Δ t \leq h^{2} /2$ stability bound — 25× more work. The adaptive integrator detected that the high-frequency modes had already decayed and the slow first mode could be tracked with much larger steps; FTCS, in contrast, has to use the same conservative $Δ t$ regardless of solution structure.

Stiff vs non-stiff: choosing the integrator

Three categories of ODE integrator suit MOL applications:

Explicit Runge-Kutta (RK4, Dormand-Prince RK5(4) = scipy's RK45). Cheap per step. Good for HYPERBOLIC-dominated problems where the stiffness is moderate and the CFL limit is the binding constraint anyway. Adaptive step control is highly effective.
Implicit methods for stiff problems: BDF (backward differentiation formulas, scipy's BDF), Rosenbrock methods, implicit Runge-Kutta. Allow $Δ t$ orders of magnitude larger than explicit methods on stiff systems, at the cost of solving a linear (or nonlinear) system per step. Essential for PARABOLIC problems on fine meshes.
Geometric / symplectic integrators: for Hamiltonian PDEs (wave equation, Schrödinger, KdV) where you want long-time conservation of energy or other invariants. Störmer-Verlet, Forest-Ruth, and higher-order symplectic schemes preserve the symplectic structure of the underlying flow. See time-dependent Schrödinger for the quantum-dynamics version.

Generalizations

MOL works for any time-dependent PDE — wave, Schrödinger, reaction-diffusion, advection-diffusion-reaction — and any spatial discretization. Three notable variants:

MOL on unstructured meshes: combine with finite-element spatial discretization. The resulting ODE system is $M u^{'} (t) = K u (t)$ with the mass matrix $M$ ; multiply through by $M^{- 1}$ (or factor it once) and integrate. This is exactly how the time-dependent FEM heat example works.
Lines in 2D and 3D: discretize transversely (e.g. in $y$ ) and treat $x$ as the evolution variable. Useful for parabolic problems with one preferred direction.
Pseudo-time continuation for steady-state problems: turn $\hat{L} u = f$ into $u_{t} + \hat{L} u = f$ and march to steady state. This converts an elliptic problem into a parabolic one; sometimes useful for nonlinear elliptic problems where Newton iteration struggles.

Finite difference methods — the most-common spatial discretization in MOL.
Spectral methods — high-accuracy spatial discretization; pairs well with high-order time integrators for smooth problems.
Differential equations — the ODE machinery (Runge-Kutta, BDF, etc.) that MOL exploits.
Operator splitting — the alternative organizing principle when a PDE has multiple distinct physical effects.
FEM heat 1D — a worked example using MOL with FE spatial discretization.