Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 10

Boundary Value Problems for Ordinary Differential Equations

The exercises below practice solver choices and numerical checks. Eight physics projects introduce new boundary value and eigenvalue models. Each includes an exact result or a physical identity for validation. Some projects require extending the supplied driver or using a library solver; those changes are stated explicitly.

Chapter 10 codes Chapter 10 resources Chapter 10 updates & errata Chapter 10 overview

Additional exercises

Use the supplied codes or a documented library routine. State whether n counts intervals or points, separate solver tolerance from grid spacing, and compare solutions on common coordinates.

  1. A manufactured solution with a first-derivative term

    Use Code 1 for y″ + 2y′ − y = 2ex on [0,1], with y(0) = 1 and y(1) = e. The exact solution is y = ex. Check the signs of both off-diagonal coefficients and the endpoint contributions.

    Use 20, 40, and 80 intervals. Report maximum and RMS errors, with the RMS definition stated explicitly. Compare with shooting after writing y″ = 2ex − 2y′ + y. Verify second-order finite-difference convergence.

  2. A boundary layer that a uniform grid may not resolve

    Solve εy″ + y′ = 0 on [0,1], with y(0) = 0 and y(1) = 1. Use ε = 0.1, 0.02, and 0.005. The exact solution is y = (1 − e−x/ε)/(1 − e−1/ε).

    Use Code 1 with p = 1/ε, q = r = 0 and 20, 100, and 400 intervals. Inspect overshoots and nonmonotone values when h/(2ε) exceeds one. Compare with a library BVP solver with a refined mesh near x = 0. Explain why a linear-system solve can succeed although the differential equation is poorly resolved.

  3. Pure derivative conditions: compatibility and a missing reference value

    Extend Code 2 to y″ = r(x) with derivatives prescribed at both ends. Its existing derivative options prescribe a derivative at one end and a value at the other; simultaneous derivative conditions require a new boundary option. Integrating the equation shows that y′(b) − y′(a) must equal the integral of r.

    On [0,1], try r = 0 with both derivatives zero, then with derivatives zero and one. Explain nonuniqueness in the first case and incompatibility in the second. For r = 0, recover a unique solution by supplying one reference value or mean value in addition to compatible derivative data. Do not interpret a singular tridiagonal system as an ordinary near-zero numerical error.

  4. Newton convergence is different from grid convergence

    Use Code 3 on 20, 40, and 80 intervals, with Newton tolerances 10−4, 10−8, and 10−10. Keep the supplied nonlinear equation and exact comparison.

    Tabulate iteration count, final discrete residual, and maximum solution error. Tighten the tolerance at fixed h, then refine h at fixed tolerance. Explain why a residual near machine precision can coexist with a visible error relative to sin(πx). Repeat with initial guesses 0 and 2sin(πx).

  5. An eigenvalue scan can miss several roots

    In Code 4, search the uniform-string problem on [0,100] with scan steps 1, 20, and 50. Keep 200 grid intervals and a bracket-width tolerance of 10−10.

    Use the known roots π2, 4π2, and 9π2 to inspect each scan interval. Record missing roots and any solver message reporting that fewer modes were found than requested; do not assume every bracket contains one root. Count interior nodes and compare with a fine scan. Explain why tightening bisection tolerance cannot recover roots never bracketed.

  6. Eigenvector signs and numerical checks

    Use Code 5 and compare its selected eigenpairs with a symmetric library eigensolver on the same matrix. Align each eigenvector’s overall sign before computing pointwise differences.

    Check the weighted overlap hΣψjψk, the residual Hψ − Eψ, and node counts. Multiply one mode by −1 and show which checks change. Normalize one mode incorrectly by omitting h and show why Euclidean vector normalization is different from wavefunction normalization.

Additional physics projects

These projects provide practice with different physical models and numerical checks. Parameter values are reproducible starting points, not measured data. For each project, report the model, discretization, convergence, numerical checks, and physical interpretation.

1. A cooling fin with a convecting tip

A uniform fin loses heat along its sides. For excess temperature θ = T − T∞, the steady equation is θ″ = m2θ, with θ(0) = 1 and θ′(1) + Bθ(1) = 0. Use dimensionless m = 2 and B = 0.5. This includes distributed cooling as well as heat loss at the tip.

Use Code 2 with the right Robin option (key 4), p = 0, q = −m2, r = 0, α = B, β = 1, and right_value = 0. Refine the grid, then vary m and B. Plot temperature and the dimensionless axial heat flow Q = −θ′.

Validation. θ(x) = [cosh(m(1 − x)) + (B/m)sinh(m(1 − x))]/[cosh m + (B/m)sinh m]. Check θ(1) ≈ 0.21418. Verify Q(0) − Q(1) = m2∫01θ dx and Q(1) = Bθ(1).

2. Electrostatics in a coaxial capacitor

In a charge-free cylindrical region, the potential obeys V″(r) + V′(r)/r = 0. Take inner radius a = 1, outer radius b = 3, V(a) = 1, V(b) = 0, and reduced permittivity ε = 1.

Use Code 1’s independent variable x for the radius. Set p(x) = 1/x and q = 0, and set the source term r(x) = 0; this r(x) denotes the source function, not the radius. Compute Er = −dV/dr. Plot V and Er and vary the radius ratio. Compare numerical charge per unit length at both cylindrical surfaces, taking their outward normals into account.

Validation. V(r) = ln(b/r)/ln(b/a), Er = 1/[r ln(b/a)], and capacitance per unit length C′ = 2π/ln(b/a) ≈ 5.71920. Check that 2πrEr is constant and that field energy per unit length equals C′/2 for unit potential difference.

3. A hanging cable with equal-height supports

For a cable with uniform weight per unit arc length w and constant horizontal tension H, let c = H/w. Its shape obeys y″ = √(1 + y′2)/c. Set supports at x = ±0.5 with y = 0, take upward y positive, and start with c = 0.5.

Adapt Code 3 with f = √(1+v2)/c (v = y′), df/dy = 0, and df/dv = v/[c√(1+v2)]. Use different guesses and values of c. Calculate the midpoint sag, total cable length, and support slopes. As an extension, determine c for a prescribed cable length by an outer scalar root search.

Validation. y(x) = c[cosh(x/c) − cosh(0.5/c)]. Sag = c[cosh(0.5/c) − 1] ≈ 0.271540 for c = 0.5, and length = 2c sinh(0.5/c) ≈ 1.175201. Verify symmetry and compare all three derived quantities with these formulas.

4. Thermal ignition and two steady temperature profiles

A dimensionless slab model of temperature-sensitive heat generation gives θ″ + δeθ = 0 on [0,1], with θ(0) = θ(1) = 0. Begin with δ = 1. The model is also called the Bratu problem.

Adapt Code 3 using f = −δey, df/dy = −δey, df/dv = 0 (v = y′). Try initial guesses 0 and 5sin(πx). Use a library solver or damped Newton steps if full Newton fails. Continue each solution branch as δ increases; distinguish convergence failure from loss of a steady solution. A steady BVP calculation alone does not establish dynamical stability.

Validation. Parameterize exact profiles by k > 0: θ(x) = 2ln[cosh(k/2)/cosh(k(x − 1/2))], with δ = 2k2/cosh2(k/2). For δ = 1 the midpoint temperatures are approximately 0.14054 and 4.09147. The branches meet near δ = 3.51383, determined by k tanh(k/2) = 2.

5. Deflection of a simply supported beam

For a uniform Euler–Bernoulli beam with downward deflection y and uniform load q, EI y(4) = q. Use L = EI = q = 1, with y(0) = y(L) = 0 and y″(0) = y″(L) = 0.

Introduce the curvature κ = y″. First solve κ″ = q/(EI), with κ(0) = κ(L) = 0; then solve y″ = κ, with zero endpoint deflections. Reuse the finite-difference construction in Code 1, feeding the first solution into the second system. Refine the grid, examine curvature, and repeat for a smooth localized load.

Validation. For uniform load, κ(x) = qx(x − L)/(2EI) and y(x) = qx(L3 − 2Lx2 + x3)/(24EI). The maximum deflection is 5qL4/(384EI) = 0.0130208333. The bending moment is proportional to EIκ; its sign depends on the chosen mechanical convention.

6. Quantum states above a reflecting surface

The dimensionless quantum bouncer has V(x) = x for x ≥ 0, ψ(0) = 0, and decay as x → ∞. Its Hamiltonian is −ψ″/2 + xψ = Eψ.

Adapt Code 5 to [0,L] with L = 12, using a library tridiagonal eigensolver for fine grids. Compute the lowest three states, refine h, then increase L at approximately fixed h. Compare the classical turning point x = E with each probability density and count nodes.

Validation. The exact energies are Ej = −aj/21/3, where aj are the negative zeros of Ai: −2.3381074105, −4.0879494441, −5.5205598281. Reference energies are approximately 1.855757, 3.244608, and 4.381671. The wavefunction is proportional to Ai(21/3(x − E)).

7. Vibrational bound states in a Morse potential

Use V(x) = D(e−2ax − 2e−ax), with D = 8, a = 1, and ℏ = mass = 1. The minimum is −D and the dissociation threshold is zero.

Adapt Code 5 on [−3,12], using a tridiagonal library solver for refined grids. Identify negative-energy localized states, refine the grid, and move the right boundary outward. Plot the first few wavefunctions and examine the unequal energy spacing. Check whether the highest bound state is more sensitive to box size.

Validation. Let s = √(2D)/a. Exact energies are En = −(a2/2)(s − n − 1/2)2, with integer n ≥ 0 and n < s − 1/2. These parameters give four bound states: −6.125, −3.125, −1.125, −0.125. Positive finite-box eigenvalues are not additional bound states.

8. A quartic quantum oscillator: scaling and the virial theorem

Consider H = −(1/2)d2/dx2 + gx4, with g = 0.25, 1, and 4. This has no harmonic-oscillator energy formula.

Adapt Code 5 and use symmetric intervals large enough for decaying tails. Compute the first three energies and normalized wavefunctions. After the change of variable x = g−1/6ξ, use the same ξ interval and grid for each g, mapping it to x for the calculation. Use the same discrete kinetic operator to evaluate kinetic and potential expectation values.

Validation. Scaling gives En(g) = g1/3En(1), and ψg(x) = g1/12ψ1(g1/6x). In the continuum, 2⟨T⟩ = 4⟨V⟩, so ⟨T⟩ = 2E/3 and ⟨V⟩ = E/3. Check that the virial defect decreases with grid and interval refinement, and that states alternate parity.

← Back to Chapter 10