Chapter 10
Boundary Value Problems for Ordinary Differential Equations
The companion programs make shooting, finite differences, and eigenvalue calculations visible step by step. The resources below provide tested library alternatives, analytic solutions and numerical benchmarks for checking results, and guidance on choosing a method.
Software and library resources
Library BVP solvers accept a first-order system, boundary residuals, an initial mesh, and a solution guess. MATLAB bvp4c and SciPy solve_bvp use fourth-order collocation with residual control; MATLAB bvp5c uses fifth-order collocation and controls an estimate of the solution error. Check each solver’s tolerance definition and reported status, then repeat with tighter tolerances and compare solutions.
On a narrow screen, scroll horizontally to view the table.
| Numerical task | MATLAB | Python | C++ | Fortran |
|---|---|---|---|---|
| Two-point BVP | bvp4c; bvp5c | SciPy solve_bvp | Use the supplied shooting or finite-difference programs, or call a Fortran BVP package through an appropriate language interface. | BVP_SOLVER / BVP_SOLVER-2 (Fortran 90/95); COLNEW and COLSYS (collocation) |
| Tridiagonal linear system | Build a sparse tridiagonal matrix with spdiags, then solve with A\rhs | solve_banded or spsolve | Call LAPACK DGTSV through a suitable C++ interface; it uses partial pivoting. | LAPACK DGTSV uses partial pivoting. |
| Symmetric matrix eigenproblem | eig; eigs for selected eigenpairs | NumPy eigh; SciPy eigh_tridiagonal | LAPACK DSTEV or DSTEVR through a suitable C++ interface. | LAPACK DSTEV for all eigenpairs; DSTEBZ for selected eigenvalues; DSTEVR for selected eigenpairs. |
| Shooting root refinement | fzero | brentq for a bracketed root; newton uses secant iteration when no derivative is supplied. Verify the endpoint residual. | Use the supplied secant or bisection routine. A bracketed search needs a valid sign change. | The supplied secant and bisection routines keep the search explicit. A bracketed method needs a valid sign change. |
Benchmark equations and program defaults
These checks use the supplied program defaults. Here n counts grid intervals, so n+1 points include the endpoints; shooting uses an adaptive grid. Iteration and step counts can vary slightly with floating-point arithmetic. Compare the solution or eigenvalues as well as the reported residual.
On a narrow screen, scroll horizontally to view the table.
| Equation and boundaries | Analytic comparison | Defaults and checks |
|---|---|---|
| 0: shootingy″ = y; y(0) = y(1) = 1 | y = cosh(x − 0.5)/cosh(0.5) | [0,1]; initial slopes 0, 1; RKF45 tolerance 10−10; boundary tolerance 10−8; hmax = 0.05. Slope ≈ −0.4621171572; 3 trials, final trial 29 accepted and 1 rejected step. |
| 1: linear finite differencesy″ − y = 0; y(0) = y(1) = 1 | Same solution as Code 0 | 20 intervals, h = 0.05; maximum grid error ≈ 2.1338 × 10−5. |
| 2: boundary choicesy″ = −sin(πx); y(0) = 0, y′(1) = 0 | y = sin(πx)/π2 + x/π | Default key 2; 20 intervals; maximum grid error ≈ 6.5477 × 10−4. Exact solutions for other keys are computed from their boundary data. |
| 3: nonlinear finite differencesy″ = y2 − sin2(πx) − π2sin(πx); y(0) = y(1) = 0 | y = sin(πx) | 20 intervals; zero initial guess; discrete-residual tolerance 10−10. 4 Newton iterations; maximum grid error ≈ 1.7496 × 10−3. |
| 4: Numerov modesy″ = −λy; y(0) = y(1) = 0 | λm = m2π2; ym = √2 sin(mπx) | 200 intervals; 3 modes; search [0,100], scan step 1, bracket width 10−10. Numerical λ: 9.8696043986, 39.4784174441, 88.8264377845. |
| 5: quantum matrix−ψ″/2 + x2ψ/2 = Eψ; ψ(±6) = 0 | Infinite-line energies Ej = j + 1/2 | 120 intervals; 119 interior unknowns; 3 states. Numerical E: 0.49968730, 1.49843574, 2.49593063. At L = 6, the error in these states is dominated by grid spacing. Check both grid and interval size for another potential. |
Choosing a boundary value method
On a narrow screen, scroll horizontally to view the table.
| Problem | Good method to try | Notes |
|---|---|---|
| Smooth nonstiff BVP | Shooting | Shooting is easy to adapt; use a library collocation solver if trial integrations amplify errors. Check the entire trajectory. |
| Unstable or rapidly growing shooting solutions | Finite differences, collocation, or multiple shooting | Use finite differences, a library collocation solver, or multiple shooting. Tightening an IVP tolerance alone may not overcome amplification of trial errors. |
| Infinite domain | Finite differences or collocation on a truncated interval | Truncate to a finite interval and check dependence on its size while keeping h approximately fixed. If the far-field decay rate is known, an asymptotic Robin condition may improve the boundary approximation. For Schrödinger eigenproblems the decay rate generally depends on the unknown energy E. |
| Linear BVP | Centered finite differences | Centered finite differences give a tridiagonal system. Resolve layers with enough points and use pivoting when needed. |
| Nonlinear BVP | Finite differences with Newton iteration | Use Newton iteration with correct partial derivatives. Try different initial guesses and continuation. Failure does not prove nonexistence. |
| Eigenvalue shooting | Numerov shooting | Avoid the trivial zero solution by fixing a nonzero amplitude. Scan conservatively, check nodes, and refine brackets. |
| Many eigenpairs | Sparse or tridiagonal matrix eigensolver | Use a matrix formulation. Preserve symmetry, omit known zero endpoints from the unknowns, and exploit sparse or tridiagonal storage. |
For smooth, resolved solutions, halving h should reduce the finite-difference errors in Codes 1–3 by about four and the uniform-string Numerov eigenvalue errors in Code 4 by about sixteen. A small iteration residual does not remove grid error. For truncated infinite domains, vary the interval size while keeping h approximately fixed. Check eigenfunction normalization and orthogonality, and use a tridiagonal eigensolver for large grids.
Selected references and documentation
- U. M. Ascher, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, SIAM, 1995.
- J. Kierzenka and L. F. Shampine, “A BVP Solver Based on Residual Control and the MATLAB PSE,” ACM Transactions on Mathematical Software 27, 299–316 (2001).
- L. F. Shampine, P. H. Muir, and H. Xu, “A User-Friendly Fortran BVP Solver,” Journal of Numerical Analysis, Industrial and Applied Mathematics 1, 201–217 (2006).
- W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007. Chapter 18 covers two-point boundary value problems, shooting, and relaxation.
- J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007. See the quantum-mechanics material for numerical Schrödinger calculations.