28 Numerical Methods for ODEs
An ordinary differential equation prescribes the derivative of an unknown function. A time stepping method replaces the exact flow by a sequence of computable updates. The numerical questions are accuracy, stability, stiffness, and cost per step.
28.1 Initial Value Problems
An initial value problem (IVP) has the form \[ \begin{align} \by'(t)=\boldsymbol{f}(t,\by(t)), \qquad \by(t_0)=\by_0, \end{align} \] where \(\by(t)\in\fR^m\). The scalar case is \(m=1\).
(IVP vs. BVP) In an IVP, all data are given at one time. In a boundary value problem, conditions are imposed at different points, such as \(u(0)=u(1)=0\). IVPs are naturally marched forward; BVPs usually lead to global linear or nonlinear systems.
For the IVP in the result above, the exact flow map \(\varphi_h\) satisfies \[ \begin{align} \by(t+h)=\varphi_h(t,\by(t)). \end{align} \] A one-step method replaces \(\varphi_h\) by a computable map \(\Phi_h\): \[ \begin{align} \by_{n+1}=\Phi_h(t_n,\by_n), \qquad t_n=t_0+nh. \end{align} \]
An ODE is autonomous if \(\boldsymbol{f}\) does not depend explicitly on time: \[ \begin{align} \by'=\boldsymbol{f}(\by). \end{align} \] Non-autonomous systems can be made autonomous by adding \(t'=1\) as one more equation.
If \(\boldsymbol{f}(t,\by)\) is continuous in \(t\) and Lipschitz continuous in \(\by\) on a neighborhood of \((t_0,\by_0)\), then the IVP in the result above has a unique local solution.
(Role of Lipschitz continuity) Numerical convergence is meaningful only when the exact IVP is well posed. If nearby initial data can lead to different solution branches, small rounding or truncation errors may not remain small.
Show that \(y'=2\sqrt{|y|}\), \(y(0)=0\), has more than one solution.
Identify which hypothesis of the result above fails.
Explain why a time stepping method cannot be expected to select a unique physical branch without extra information.
28.2 Error and Convergence
For a one-step method \(\by_{n+1}=\Phi_h(t_n,\by_n)\), the local truncation error is the one-step defect made from exact data: \[ \begin{align} \boldsymbol{\tau}_{n+1} = \frac{\by(t_{n+1})-\Phi_h(t_n,\by(t_n))}{h}. \end{align} \] The method is consistent if \(\|\boldsymbol{\tau}_{n+1}\|\to 0\) as \(h\to 0\).
The global error at \(t_n\) is \[ \begin{align} \be_n=\by(t_n)-\by_n. \end{align} \] A method has order \(p\) on a fixed interval \([t_0,T]\) if \[ \begin{align} \max_{0\leq n\leq N}\|\be_n\|=O(h^p), \qquad Nh=T-t_0. \end{align} \]
For a stable one-step method, a one-step defect \(O(h^{p+1})\) gives global error \(O(h^p)\) on a fixed time interval.
Suppose each step contributes a local error of size \(O(h^{p+1})\).
Count the number of steps needed to reach a fixed final time.
Explain why this gives \(O(h^p)\) global error when previous errors are not amplified strongly.
State where stability enters this argument.
For time stepping methods, consistency means the discrete update matches the differential equation as \(h\to 0\). Stability means perturbations introduced at one step do not grow uncontrollably over a fixed interval. Convergence means \(\by_n\to\by(t_n)\) as \(h\to 0\).
For a consistent one-step method applied to a Lipschitz IVP, stability on a fixed interval implies convergence. For linear multistep methods, consistency plus zero-stability implies convergence.
Write the error recurrence for a one-step method by subtracting the numerical update from the exact one-step update.
Bound the propagated error using a Lipschitz constant for \(\boldsymbol{f}\).
Add the local defect term from the result above.
Apply a discrete Gronwall argument to justify the first statement in the result above.
28.3 Euler Methods
The explicit Euler method uses the slope at the left endpoint: \[ \begin{align} \by_{n+1}=\by_n+h\boldsymbol{f}(t_n,\by_n). \end{align} \] It is explicit because the right side is known before step \(n+1\) is computed.
Forward Euler has local truncation error \(O(h)\) in the sense of the result above, one-step defect \(O(h^2)\), and global error \(O(h)\) for smooth solutions on fixed intervals.
Taylor expand the exact solution: \[\by(t_{n+1})=\by(t_n)+h\by'(t_n)+\frac{h^2}{2}\by''(\xi_n).\]
Substitute \(\by'(t_n)=\boldsymbol{f}(t_n,\by(t_n))\).
Compare with the result above.
Use the result above to obtain the global order.
(Forward Euler on decay) For \(y'=-y\), \(y(0)=1\), Forward Euler gives \[ \begin{align} y_n=(1-h)^n. \end{align} \] At \(t=1\) with \(h=1/4\), this gives \((3/4)^4=81/256\approx 0.3164\), while the exact value is \(e^{-1}\approx 0.3679\).
The implicit Euler method uses the slope at the right endpoint: \[ \begin{align} \by_{n+1}=\by_n+h\boldsymbol{f}(t_{n+1},\by_{n+1}). \end{align} \] At each step one solves a nonlinear equation for \(\by_{n+1}\).
For Backward Euler, define \[ \begin{align} \bg(\bY)=\bY-\by_n-h\boldsymbol{f}(t_{n+1},\bY). \end{align} \] The next step is the root \(\bg(\bY)=\bzero\). Newton’s method requires linear systems with \[ \begin{align} \bI-h\bJ_f(t_{n+1},\bY), \end{align} \] where \(\bJ_f\) is the Jacobian of \(\boldsymbol{f}\) with respect to \(\by\).
Backward Euler is first order: its one-step defect is \(O(h^2)\) and its global error is \(O(h)\) for smooth solutions.
Taylor expand \(\by(t_n)\) about \(t_{n+1}\).
Rearrange the expansion to express \(\by(t_{n+1})\) in terms of \(\by(t_n)\) and \(\by'(t_{n+1})\).
Substitute \(\by'(t_{n+1})=\boldsymbol{f}(t_{n+1},\by(t_{n+1}))\).
Compare with the result above.
The implicit trapezoidal method averages endpoint slopes: \[ \begin{align} \by_{n+1} = \by_n+\frac{h}{2}\left[\boldsymbol{f}(t_n,\by_n)+\boldsymbol{f}(t_{n+1},\by_{n+1})\right]. \end{align} \] It is also called Crank-Nicolson when used for time-dependent PDE discretizations.
The trapezoidal method has one-step defect \(O(h^3)\) and global error \(O(h^2)\) for smooth solutions.
Write the exact integral identity \[\by(t_{n+1})=\by(t_n)+\int_{t_n}^{t_{n+1}}\boldsymbol{f}(t,\by(t))\,dt.\]
Approximate the integral by the scalar trapezoidal rule from the result above.
Use the result above to explain the one-step defect.
Convert the one-step defect to global order using the result above.
28.4 Absolute Stability and Stiffness
The scalar test equation is \[ \begin{align} y'=\lambda y, \qquad \operatorname{Re}(\lambda)<0. \end{align} \] For a one-step method applied to this equation, the update has the form \[ \begin{align} y_{n+1}=R(z)y_n, \qquad z=h\lambda. \end{align} \] The function \(R\) is the stability function.
The absolute stability region of a one-step method is \[ \begin{align} \mathcal{S}=\{z\in\fC:\ |R(z)|\leq 1\}. \end{align} \] The method is stable on the test equation when \(h\lambda\in\mathcal{S}\).
For the test equation in the result above, \[ \begin{align} R_{\mathrm{FE}}(z)=1+z,\qquad R_{\mathrm{BE}}(z)=\frac{1}{1-z},\qquad R_{\mathrm{TR}}(z)=\frac{1+z/2}{1-z/2}. \end{align} \] Forward Euler is stable in the disk \(|1+z|\leq 1\). Backward Euler and trapezoidal rule are A-stable.
For Forward Euler, \(y_{n+1}=y_n+h\lambda y_n=(1+z)y_n\). For Backward Euler, \[ \begin{align} y_{n+1}=y_n+h\lambda y_{n+1} \quad\Rightarrow\quad y_{n+1}=\frac{1}{1-z}y_n. \end{align} \] For trapezoidal rule, \[ \begin{align} y_{n+1}=y_n+\frac{z}{2}(y_n+y_{n+1}) \quad\Rightarrow\quad y_{n+1}=\frac{1+z/2}{1-z/2}y_n. \end{align} \] If \(\operatorname{Re}(z)\leq 0\), then \(|1|\leq |1-z|\) and \(|1+z/2|\leq |1-z/2|\), so Backward Euler and trapezoidal rule are A-stable.
Draw the Forward Euler stability disk from \(|1+z|\leq 1\).
For real \(\lambda<0\), show that Forward Euler requires \(h\leq 2/|\lambda|\).
Verify \(|1+z/2|\leq |1-z/2|\) when \(\operatorname{Re}(z)\leq 0\) by squaring both sides.
Explain why A-stability does not imply L-stability.
A method is A-stable if its stability region contains the full left half-plane: \[ \begin{align} \{z\in\fC:\operatorname{Re}(z)\leq 0\}\subseteq\mathcal{S}. \end{align} \] An A-stable method is L-stable if additionally \(R(z)\to 0\) as \(z\to -\infty\) along the real axis.
(Damping) Backward Euler is L-stable, so very fast decaying modes are strongly damped. The trapezoidal method is A-stable but not L-stable, since \(R_{\mathrm{TR}}(z)\to -1\) as \(z\to -\infty\).
An IVP is stiff when stability forces an explicit method to take much smaller steps than accuracy requires. For a linear system \[ \begin{align} \by'=\bA\by, \end{align} \] the eigenvalues of \(\bA\) determine the explicit stability restriction through \(h\lambda_i\in\mathcal{S}\).
(Fast decay with slow forcing) The equation \[ \begin{align} y'=-1000(y-\cos t)-\sin t \end{align} \] has exact solution \(y(t)=\cos t\) for consistent initial data. The solution varies on an \(O(1)\) time scale, but Forward Euler is restricted by the fast stable mode near \(\lambda=-1000\) and requires \(h\lesssim 0.002\).
Verify directly that \(y(t)=\cos t\) solves the stiff example above.
Linearize the equation around the exact solution and identify the error equation.
Apply the result above to obtain the Forward Euler stability restriction.
Explain why this restriction is unrelated to the visible time scale of \(\cos t\).
28.5 Runge-Kutta Methods
An \(s\)-stage explicit Runge-Kutta method has stages \[ \begin{align} \boldsymbol{k}_i = \boldsymbol{f}\left(t_n+c_i h,\, \by_n+h\sum_{j<i} a_{ij}\boldsymbol{k}_j\right), \qquad i=1,...,s, \end{align} \] and update \[ \begin{align} \by_{n+1}=\by_n+h\sum_{i=1}^s b_i\boldsymbol{k}_i. \end{align} \] The coefficients are collected in a Butcher tableau.
The two-stage explicit midpoint method is \[ \begin{align} \boldsymbol{k}_1 &= \boldsymbol{f}(t_n,\by_n),\\ \boldsymbol{k}_2 &= \boldsymbol{f}(t_n+h/2,\by_n+h\boldsymbol{k}_1/2),\\ \by_{n+1} &= \by_n+h\boldsymbol{k}_2. \end{align} \] It has global order \(2\).
Taylor expand the exact solution to second order.
Taylor expand \(\boldsymbol{f}(t_n+h/2,\by_n+h\boldsymbol{k}_1/2)\) about \((t_n,\by_n)\).
Match terms to show second-order consistency.
Derive the stability function of explicit midpoint on \(y'=\lambda y\).
The classical fourth-order Runge-Kutta method is \[ \begin{align} \boldsymbol{k}_1 &= \boldsymbol{f}(t_n, \by_n),\\ \boldsymbol{k}_2 &= \boldsymbol{f}(t_n + h/2, \by_n + h \boldsymbol{k}_1 / 2),\\ \boldsymbol{k}_3 &= \boldsymbol{f}(t_n + h/2, \by_n + h \boldsymbol{k}_2 / 2),\\ \boldsymbol{k}_4 &= \boldsymbol{f}(t_n + h, \by_n + h \boldsymbol{k}_3),\\ \by_{n+1} &= \by_n + \frac{h}{6}(\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4). \end{align} \] It has one-step defect \(O(h^5)\) and global error \(O(h^4)\) for smooth nonstiff problems.
Applied to \(y'=\lambda y\), RK4 has stability function \[ \begin{align} R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}. \end{align} \] This is the degree-four Taylor polynomial of \(e^z\).
Apply the result above to \(y'=\lambda y\).
Express each stage in terms of \(z=h\lambda\) and \(y_n\).
Substitute the stages into the RK4 update.
Derive the stability polynomial in the result above.
(High order is not stiffness control) RK4 has higher accuracy than Forward Euler on nonstiff problems, but its stability region is still bounded. For stiff decay, an explicit high-order method can remain limited by stability rather than accuracy.
28.6 Adaptive Time Stepping
An embedded pair computes two approximations of different orders from shared stages: \[ \begin{align} \by_{n+1}^{[p]}, \qquad \by_{n+1}^{[p+1]}. \end{align} \] The difference \[ \begin{align} \be_{n+1}=\by_{n+1}^{[p+1]}-\by_{n+1}^{[p]} \end{align} \] estimates the local error of the lower-order method.
Adaptive solvers usually accept a step when \[ \begin{align} \left\| \frac{\be_{n+1}} {\mathrm{atol}+\mathrm{rtol}\,|\by_{n+1}|} \right\| \leq 1, \end{align} \] where the absolute value and division are componentwise.
If the local error estimate satisfies \(\|\be\|\approx Ch^{p+1}\), then a step size targeting tolerance \(\varepsilon\) is \[ \begin{align} h_{\mathrm{new}} = s\,h_{\mathrm{old}} \left(\frac{\varepsilon}{\|\be\|}\right)^{1/(p+1)}, \end{align} \] where \(0<s<1\) is a safety factor.
Assume \(\|\be_{\mathrm{old}}\|\approx C h_{\mathrm{old}}^{p+1}\).
Solve \(C h_{\mathrm{new}}^{p+1}\approx \varepsilon\).
Eliminate \(C\) using the old error estimate.
Explain why production codes also limit the maximum growth and shrinkage of \(h\).
(RK45) Dormand-Prince RK45 is a common embedded explicit pair. It is a good default for smooth nonstiff IVPs, not for stiff systems.
28.7 Linear Multistep Methods
A \(k\)-step linear multistep method has the form \[ \begin{align} \sum_{j=0}^k \alpha_j \by_{n+j} = h\sum_{j=0}^k \beta_j \boldsymbol{f}(t_{n+j},\by_{n+j}). \end{align} \] It reuses previous solution values and previous right-hand side evaluations.
The two-step Adams-Bashforth method is explicit: \[ \begin{align} \by_{n+1} = \by_n+\frac{h}{2}\left(3\boldsymbol{f}(t_n,\by_n)-\boldsymbol{f}(t_{n-1},\by_{n-1})\right). \end{align} \] It has global order \(2\).
The two-step Adams-Moulton method is implicit: \[ \begin{align} \by_{n+1} = \by_n+\frac{h}{2}\left(\boldsymbol{f}(t_{n+1},\by_{n+1})+\boldsymbol{f}(t_n,\by_n)\right). \end{align} \] This is the trapezoidal method from the result above.
The second-order backward differentiation formula is \[ \begin{align} \frac{3\by_{n+1}-4\by_n+\by_{n-1}}{2h} = \boldsymbol{f}(t_{n+1},\by_{n+1}). \end{align} \] BDF methods are standard implicit methods for stiff IVPs.
A linear multistep method is zero-stable if perturbations in starting values remain uniformly bounded as \(h\to 0\) on a fixed interval. Algebraically, the roots of \[ \begin{align} \rho(\xi)=\sum_{j=0}^k \alpha_j \xi^j \end{align} \] must satisfy \(|\xi|\leq 1\), and roots on the unit circle must be simple.
For a consistent linear multistep method applied to a well-posed IVP, convergence is equivalent to zero-stability.
For AB2, write the polynomial \(\rho(\xi)\).
Check the root condition in the result above.
Repeat for BDF2.
Explain why consistency alone is not enough for a multistep method.
(BDF for stiffness) BDF methods trade implicit linear or nonlinear solves for much larger stable time steps on stiff decay problems. Variable-order BDF methods are used in many stiff ODE codes.
28.8 Systems and Eigenvalues
A higher-order scalar equation can be rewritten as a first-order system. For example, \[ \begin{align} q''=F(q,q',t) \end{align} \] becomes \[ \begin{align} q'=v,\qquad v'=F(q,v,t). \end{align} \]
For \[ \begin{align} \by'=\bA\by,\qquad \by(0)=\by_0, \end{align} \] the exact solution is \[ \begin{align} \by(t)=e^{t\bA}\by_0. \end{align} \] If \(\bA=\bV\bD\bV^{-1}\), then \[ \begin{align} \by(t)=\bV e^{t\bD}\bV^{-1}\by_0. \end{align} \] Each eigenvalue controls one modal growth, decay, or oscillation rate.
Differentiate \(e^{t\bA}\by_0\) to verify the first formula in the result above.
Substitute \(\bA=\bV\bD\bV^{-1}\) into the matrix exponential.
For \(\lambda=a+ib\), interpret \(e^{\lambda t}\) in terms of growth and oscillation.
Explain why eigenvalues of the Jacobian predict local stiffness for nonlinear systems.
(Jacobian spectrum) Near a state \(\by_*\), nonlinear dynamics are locally approximated by \(\be'=\bJ_f(\by_*)\be\). Large negative eigenvalues produce fast stable transients. Eigenvalues near the imaginary axis produce oscillation or slow decay.
28.9 Hamiltonian and Symplectic Methods
For position \(\bq\) and velocity or momentum \(\bv\), a separable Hamiltonian has \[ \begin{align} H(\bq,\bv)=\frac{1}{2}\|\bv\|^2+V(\bq), \qquad \bq'=\bv,\qquad \bv'=\bF(\bq)=-\nabla V(\bq). \end{align} \] The exact flow conserves \(H\).
For the Hamiltonian system in the result above, Velocity Verlet is \[ \begin{align} \bv_{n+1/2} &= \bv_n+\frac{h}{2}\bF(\bq_n),\\ \bq_{n+1} &= \bq_n+h\bv_{n+1/2},\\ \bv_{n+1} &= \bv_{n+1/2}+\frac{h}{2}\bF(\bq_{n+1}). \end{align} \] It is second order, time reversible, and symplectic.
For \(q'=v\), \(v'=-q\), Forward Euler maps the energy \[ \begin{align} E(q,v)=\frac{1}{2}(q^2+v^2) \end{align} \] to \((1+h^2)E(q,v)\) after one step. Thus Forward Euler has artificial energy growth on this oscillator.
Forward Euler gives \[ \begin{align} q_{n+1}=q_n+hv_n,\qquad v_{n+1}=v_n-hq_n. \end{align} \] Therefore \[ \begin{align} q_{n+1}^2+v_{n+1}^2 = (q_n+hv_n)^2+(v_n-hq_n)^2 = (1+h^2)(q_n^2+v_n^2). \end{align} \] Multiplying by \(1/2\) proves the claim.
Apply Velocity Verlet to \(q'=v\), \(v'=-q\).
Show that reversing the step size from \(h\) to \(-h\) undoes one Verlet step.
Compare long-time energy behavior of Forward Euler, RK4, and Velocity Verlet numerically.
(Symplectic accuracy) For long Hamiltonian simulations, preserving phase-space geometry can matter more than minimizing one-step error. Symplectic methods usually keep the energy error bounded and oscillatory over long times.
28.10 Method Choice
(SciPy interface) In scipy.integrate.solve\_ivp, methods such as RK45 and DOP853 are explicit nonstiff solvers. Methods such as BDF and Radau are intended for stiff problems and benefit from Jacobian information.
For \(y'=-y\), compare Forward Euler, RK4, and RK45 on accuracy versus function evaluations.
For the stiff example in the result above, compare RK45 and BDF using the same tolerance.
For a sparse semidiscrete heat equation \(\by'=\bA\by\), identify why Jacobian sparsity matters.
For the harmonic oscillator, compare final-time error and energy error for RK4 and Velocity Verlet.