Method of Lines and Time Integration

The method of lines freezes the spatial discretization and treats the remaining time dependence as a large ODE system. This modularity is powerful only when spatial and temporal stability are compatible. A high-order Runge–Kutta method cannot rescue an unstable spatial operator, and a stable spatial discretization can be destabilized by choosing a timestep outside the integrator’s stability region.

Required background. Ordinary Differential Equations supplies Runge–Kutta, multistep, implicit, and adaptive time integration.

Helpful background. Formulations of General Relativity supplies the continuum evolution systems and constraints, while Partial Differential Equations and Green Functions provides characteristics, spatial operators, and boundary conditions.

From a PDE to a semidiscrete ODE

After choosing a grid or element space, write

$$\frac{dU}{dt}=L_h(U,t), \qquad U(0)=U_h^0.$$

The operator $L_h$ includes spatial derivatives, numerical fluxes, sources, dissipation, and numerical boundary conditions. The time integrator sees this complete operator, not the continuum differential equation. For linear $U'=A_hU$, stability depends on the eigenvalues and, for nonnormal $A_h$, may also depend on transient resolvent growth.

An explicit Runge–Kutta method advances $U^{n+1}=R(\Delta t A_h)U^n$. Plotting or tabulating the semidiscrete spectrum against the region $\lvert R(z)\rvert\le1$ is a direct compatibility test. The derivation and preservation properties of strong-stability-preserving Runge–Kutta methods begin in Gottlieb, Shu, and Tadmor 2001, § 1, pp. 89–90, and § 2.1, Eqs. 2.1–2.8.

A decisive unstable pairing

For periodic advection $u_t+u_x=0$ on $0\le x\lt 2\pi$, discretize with $(D_0u)_j=(u_{j+1}-u_{j-1})/(2h)$. Its semidiscrete eigenvalues are imaginary. Forward Euler has $R(z)=1+z$, so for any nonzero real $\alpha$,

$$\lvert R(i\alpha)\rvert^2=1+\alpha^2\gt 1.$$

Thus every mode with a nonzero derivative eigenvalue is amplified by Forward Euler. At fixed nonzero Courant number $\Delta t/h$, the scheme is unstable as the grid is refined: high-frequency amplification over a fixed physical time has no resolution-independent bound. A much more restrictive scaling $\Delta t=O(h^2)$ can bound growth on a fixed interval, but does not make Forward Euler contractive on imaginary eigenvalues and discards the usual hyperbolic timestep advantage. As a negative test, evolve $u(0,x)=\sin x+10^{-12}\sin(Nx/4)$ on $N=64,128,256$ grids with $\Delta t/h=0.1$ to $t=20\pi$ and record modal amplification. This detects a harness that reports accuracy only at very short final times.

The positive benchmark replaces Forward Euler by classical fourth-order Runge–Kutta. On the same domain use $N=64,128,256$, $\Delta t/h=1$, and final time $t=2\pi$, shortening only the last step to land exactly on the endpoint. Evolve $u(0,x)=\sin x$, compare with $u(t,x)=\sin(x-t)$, and report phase, amplitude, $L^2$, and $L^\infty$ errors. The expected production order is two because $D_0$ limits the fourth-order time integrator.

Separating spatial and temporal order

For a temporal test, compare to the exact semidiscrete solution rather than the continuum wave. Since $D_0e^{ix_j}=i(\sin h/h)e^{ix_j}$, the nodal ODE has the exact solution

$$u_j^{\rm semi}(t)=\sin\!\left(x_j-\frac{\sin h}{h}t\right).$$

At fixed $N=64$, use $\Delta t/h=2,1,1/2,1/4$ and final time $2\pi$. The modal error relative to this solution decreases by approximately $16$ each time the step is halved, exposing RK4’s fourth order without a spatial-error floor. Equivalently, compare $R(-i\Delta t\sin h/h)^{n}$ with $e^{-it\sin h/h}$ when $n\Delta t=t$. In contrast, comparing a very fine-grid RK4 run directly to $\sin(x-t)$ can already be dominated by the $O(h^2)$ spatial phase error at every tested timestep.

For spatial order, compare against the continuum wave on $N=64,128,256$ with $\Delta t/h=1/16$; for the combined scheme use the same grids at fixed $\Delta t/h=1$. These three comparisons distinguish the time integrator, the spatial derivative, and their production pairing.

For shock-capturing operators, strong-stability-preserving methods preserve a convex stability property enjoyed by Forward Euler under a reduced timestep. The relevant decomposition and timestep coefficient are stated in Gottlieb, Shu, and Tadmor 2001, § 2.2, Eq. 2.9, Lemma 2.1, and Eq. 2.10. “SSP” does not mean unconditionally stable; it transfers a forward-Euler property under specified assumptions and a CFL bound.

Explicit, implicit, and IMEX choices

Explicit Runge–Kutta methods are simple and parallel but limited by the fastest resolved wave. Implicit methods enlarge the stability region at the cost of nonlinear and linear solves. An implicit–explicit method treats stiff sources implicitly and transport explicitly, but the splitting and stage coupling determine its order and stability. Additive order conditions and distinct IMEX tableaux are given in Kennedy and Carpenter 2003, §§ 2–4 and Tables 3–5. Pin the exact tableau; a method name such as “third-order IMEX” is not unique.

For relativistic radiation or resistive systems, test the zero-stiffness limit, the stiff relaxation limit, and an intermediate regime. A method can pass smooth advection while losing its designed order as the stiffness parameter changes. Adaptive local error control must also respect global CFL and synchronization constraints on refined meshes.

Environment, rerun, and failure criteria

Pin the implementation revision, spatial operator, boundary closure, dissipation, integrator tableau, precision, timestep controller, Courant prescription, stage boundary updates, and norm definitions. A clean rerun exports the semidiscrete spectrum or a justified bound, stability-function comparison, exact-wave errors, spatial and temporal refinement tables, constraint histories, and work per stage.

The validity regime assumes the spatial spectrum or energy estimate used in selecting $\Delta t$ represents the active mesh and boundaries. Failure modes include Forward Euler growth with a centered operator, temporal order reduction, stage-inconsistent boundary data, a CFL violation masked by dissipation, and a stiff source handled explicitly. Reject convergence measured only at fixed timestep, because spatial refinement then changes the Courant number and confounds the comparison.

Check your understanding

A fourth-order spatial operator combined with second-order Runge–Kutta converges at second order when grid and timestep are halved together. Which component limits the result, and how can that diagnosis be verified?

Solution

The temporal integrator is the likely limiter at fixed Courant number. Hold a fine grid fixed and halve the timestep to confirm second-order temporal convergence; then use a much smaller timestep and refine space to recover the spatial order. Both sequences are needed because pre-asymptotic spatial error could otherwise mimic the same rate.

References

• Gottlieb, S., Shu, C.-W., and Tadmor, E. (2001), “Strong Stability-Preserving High-Order Time Discretization Methods,” SIAM Review 43, 89–112. DOI.
• Kennedy, C. A., and Carpenter, M. H. (2003), “Additive Runge–Kutta Schemes for Convection–Diffusion–Reaction Equations,” Applied Numerical Mathematics 44, 139–181. DOI.