Chapter 9
Initial Value Problems for Ordinary Differential Equations
The companion programs introduce fixed-step methods, adaptive RKF45, first-order systems, oscillations, projectile motion, and velocity Verlet for an orbit. The resources below extend the ODE routines discussed in Appendix B with documentation, exact solutions for checking the programs, and practical guidance on accuracy and method selection. The supplied programs are Codes 0–5; the stiff solvers listed here are for further study.
Software and library resources
A tested solver is useful when moving beyond the educational programs. Supply the derivatives, initial values, integration interval, and tolerances. Write a second-order equation as first-order equations for position and velocity. Check the solver’s return status and compare with an exact solution or a more accurate calculation.
On a narrow screen, scroll horizontally to view all language columns.
| Numerical task | MATLAB | Python | C++ | Fortran |
|---|---|---|---|---|
| Nonstiff scalar ODE or first-order systemAdaptive explicit Runge–Kutta integration. | ode45 |
solve_ivp with method='RK45' |
Boost.Odeint: runge_kutta_dopri5 with a controlled stepper |
DOPRI5; source and example driver |
| Nonstiff problem requiring high accuracyCompare cost after tightening tolerances. | ode113, a variable-order Adams solver |
solve_ivp with method='DOP853' |
Boost.Odeint: controlled runge_kutta_fehlberg78 |
DOP853; source and example driver |
| Stiff initial-value problemUse a method designed for stiffness. | ode15s |
solve_ivp with method='Radau' or 'BDF'; 'LSODA' switches between Adams and BDF |
Boost.Odeint: rosenbrock4_controller; requires Jacobian information and Boost.uBLAS state vectors |
ODEPACK: DLSODE with its stiff option, or DLSODA for automatic switching; RADAU5 is another option |
| Conservative particle dynamicsAcceleration depends only on position. | Start with Code 5: fixed-step velocity Verlet, available in all four languages. For advanced structure-preserving methods, see Hairer’s geometric integration codes. | |||
MATLAB ode45, SciPy RK45, and
DOPRI5 use Dormand–Prince formulas, whereas the
companion programs use Fehlberg’s RKF45 formulas. Similar method
names do not imply identical coefficients, error tests, accepted
steps, or results. ODEPACK supplies both single- and double-precision
routines; the prefix D denotes double precision.
Appendix B also discusses Julia. Readers using that language can
consult the SciML ODE solver documentation
for nonstiff and stiff methods, including Tsit5.
Benchmark equations and exact solutions
These solutions satisfy the equations and initial values shown below. Compare at the same times, including intermediate points. An exact solution applies to the stated model and parameters; changing the force or initial data may require a different reference solution. Rows for Codes 0–5 use the default model, initial values, and integration interval in the supplied programs.
On a narrow screen, scroll horizontally to view the full benchmark table.
| Problem and initial data | Exact solution | Why it is useful |
|---|---|---|
Code 0: x′ = x − t2 + 1, x(0) = 0.5 |
x(t) = (t + 1)2 − 0.5et |
Compare Euler, Heun, and RK4 on 0 ≤ t ≤ 4. Halve the step and compare errors at common times. |
Code 1: x′ = 100/(1 + 10000t2), x(−1) = arctan(−100); −1 ≤ t ≤ 1 |
x(t) = arctan(100t) |
A smooth but rapid transition near t = 0. Inspect how adaptive steps shrink and then grow. |
Code 2: x1′ = x2, x2′ = −x1, x3′ = −2x3, x4′ = x1 + x3; x(0) = (1,0,1,0); 0 ≤ t ≤ 10 |
x1 = cos tx2 = −sin tx3 = e−2tx4 = sin t + (1 − e−2t)/2 |
A mathematical system with oscillating, decaying, and coupled components. Check every component separately. |
Code 3: x′ = v, v′ = −x − 0.2v; x(0) = 1, v(0) = 0; 0 ≤ t ≤ 20 |
With ω = sqrt(0.99):x(t) = e−0.1t[cos(ωt) + (0.1/ω)sin(ωt)]v(t) = −e−0.1tsin(ωt)/ω |
The default underdamped spring, with m = k = 1 and damping b = 0.2. Test both position and velocity. For b = 0, the corresponding solution is x = cos t, v = −sin t. |
Code 4: no air resistance; g = 9.81, (x,y)(0) = (0,0), (vx,vy)(0) = (15,20); 0 ≤ t ≤ 4 |
x = 15ty = 20t − 4.905t2vx = 15vy = 20 − 9.81t |
Checks the ordering of four state variables and the sign of gravity. With quadratic drag, use a convergence comparison; these formulas no longer apply. |
Code 5: r″ = −r/|r|3; r(0) = (1,0), v(0) = (0,1), GM = 1; 0 ≤ t ≤ 62.83, fixed step 0.01 (6283 steps) |
r(t) = (cos t, sin t)v(t) = (−sin t, cos t)Period: 2π |
A circular orbit for checking velocity Verlet. Compare orbital phase and radius over several periods, as well as energy and angular momentum. |
Optional stiff benchmark: x′ = −1000(x − cos t) − sin t, x(0) = 0 |
x(t) = cos t − e−1000t |
A rapid initial transient followed by slow variation. Compare an explicit adaptive solver with a stiff solver at the same tolerances. This is an additional benchmark. The Chapter 9 projects page includes a related stiff-system exercise. |
Choosing a method and checking accuracy
| Problem or goal | Good starting point | What to check |
|---|---|---|
| Learn the effect of step size | Fixed-step Euler, Heun, and RK4 | For a smooth problem in the asymptotic convergence range, halving the step reduces global error by about 2, 4, and 16, respectively. Compare at a fixed final time. |
| Solve a nonstiff scalar ODE or system | Adaptive RKF45 or a tested explicit Runge–Kutta solver | The embedded estimate controls local error, not a bound on accumulated global error. Tighten tolerances and compare solutions at common times. |
| Components have different sizes or units | Scale the variables or use componentwise absolute tolerances | A single absolute tolerance is simple for teaching. Library solvers usually combine absolute and relative tolerances; choose scales appropriate to each component. |
| Explicit steps remain very small | Investigate stiffness and try a stiff solver | Small steps can also result from rapid variation, a singularity, or demanding tolerances. Compare accuracy and work before attributing the behavior to stiffness. |
| Follow a conservative orbit for many periods | Fixed-step velocity Verlet | Check phase, energy, and angular momentum. Energy is not exactly conserved. For a suitably small fixed step and a smooth, bounded orbit, its error typically remains bounded and oscillatory over long times rather than showing secular drift. Phase error can still accumulate. Halving the step should reduce global trajectory error by about 4 in the convergence range. |
| Include damping or air resistance | First-order position and velocity system with adaptive RKF45 | The standard velocity Verlet update used in Code 5 assumes acceleration depends only on position. Velocity-dependent forces require a suitable different treatment. |
Output times, plotting, and physical checks
An accurate adaptive solution can have widely spaced accepted points.
The spring and projectile companion codes limit the maximum step
with h_max to provide enough points for plotting.
A maximum step is a sampling choice as well as an integration limit;
it does not replace the error tolerance.
Library solvers can provide interpolated output at selected times:
use t_eval or dense output in SciPy, and a vector of
requested times or deval in MATLAB. These output times
need not equal the internal integration steps. Connecting a few
accepted points with straight lines is not the same as the
solver’s continuous interpolant.
The projectile companion program integrates to the requested final
time. To stop at impact, an extension can locate the downward
crossing of y = 0; MATLAB and SciPy provide event
functions for this purpose. Avoid treating the launch at
y(0) = 0 as the impact event.
For an undamped spring, check
E = mv2/2 + kx2/2.
With damping b, the physical energy satisfies
dE/dt = −bv2.
For the orbit, check specific energy
ε = |v|2/2 − GM/|r|
and specific angular momentum
ℓ = xvy − yvx.
For the default circular orbit, their exact values are
ε = −0.5 and ℓ = 1.
Conservation checks complement a trajectory comparison; a small
energy error alone does not guarantee a small phase error.
Selected references and documentation
- E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd edition, Springer, 1993. Runge–Kutta methods, error control, and multistep methods.
- E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd edition, Springer, 1996. Further study of stiffness and implicit methods.
- E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edition, Springer, 2006. Long-time integration of conservative systems.
- MATLAB: choosing an ODE solver,
odeset, anddeval - SciPy
solve_ivp: methods, tolerances, sampled output, and events - Boost.Odeint: available steppers and error control
- Hairer’s Fortran and MATLAB codes: nonstiff, stiff, and geometric integrators, with example drivers
- LLNL ODEPACK overview and Netlib ODEPACK sources. Read the solver prologue for argument definitions, workspace requirements, and return codes.
- SciML: ODE solvers for Julia, also discussed in Appendix B