Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 9

Initial Value Problems for Ordinary Differential Equations

These additional exercises and projects extend Chapter 9 without repeating its assigned problems. The exercises focus on practical questions that arise when using an ODE solver. The projects introduce different physical models, from decay chains and circuits to rotating bodies and quantum dynamics. All parameters are starting points for reproducible calculations, not measured data.

Chapter 9 codes Chapter 9 resources Chapter 9 updates & errata Chapter 9 overview

Additional exercises

Use the companion programs or a documented library solver. State the error test used in each run. Compare solutions at common times when comparing different integrations.

  1. Checking a system with a prescribed exact solution

    Solve the two-component system on 0 ≤ t ≤ 10, with u(0) = 0 and v(0) = 1:

    u′ = 2v + cos t − 2 cos(2t),
    v′ = −3u − 2 sin(2t) + 3 sin t.

    The exact solution is u = sin t, v = cos(2t). Verify it by substitution before programming the derivatives. Use Code 2 as a starting point, setting n = 2. Run with absolute tolerances 10−4, 10−7, and 10−10. Report the maximum error in each component at accepted times, final errors, and accepted and rejected steps. Explain why a small final error can conceal larger earlier errors.

  2. A change of units changes an absolute error test

    Integrate u′ = cos t and v′ = 10−5 cos(10t), starting from u(0) = v(0) = 0, to t = 10. The exact solutions are u = sin t and v = 10−6 sin(10t).

    Use Code 2 with a common absolute tolerance of 10−6, then repeat after replacing v by w = 106v, so that w′ = 10 cos(10t). Use the same initial trial step, 0.5, in both runs. Record accepted and rejected steps. Compare errors normalized by each component’s characteristic amplitude. If using a library solver, try componentwise absolute tolerances (10−8, 10−14) with relative tolerance 10−6. Explain why pointwise relative errors are unsuitable near zero crossings. Which error test preserves comparable accuracy after a change of units?

  3. An adaptive solver can miss a narrow forcing pulse

    Consider a smooth localized pulse on 0 ≤ t ≤ 1, starting from y(0) = 0:

    y′ = exp[−((t − tc)/σ)2]/(σ√π),   tc = 0.53, σ = 0.01.

    The exact solution is y(t) = [erf((t − tc)/σ) + erf(tc/σ)]/2. Use absolute tolerance 10−6 and initial step 0.1. Compare an unrestricted adaptive run with maximum steps 0.02 and 0.002; add a maximum-step limit to Code 1 if needed. Inspect the derivative as well as the solution. Repeat with tc = 0.57. Does tightening the tolerance alone reliably resolve a feature that the stage evaluations may never sample? Report what happens rather than assuming every unrestricted run fails.

  4. A known change in the differential equation

    Solve y′ = a(t)y, y(0) = 1, to t = 3, where a = 1 for t < 1 and a = −2 for t ≥ 1. The exact solution is et up to t = 1, followed by e exp[−2(t − 1)].

    Compare a single adaptive run with two integrations that end and restart at t = 1, carrying the final value into the second segment. Use absolute tolerances 10−4 and 10−8. Plot accepted step sizes near the switch. Explain why the state stays continuous although its derivative changes, and why stepping exactly to a known switch is useful.

  5. An accurate solution and an inadequate plot

    Solve y′ = 3t2, y(0) = 0, on [0,1] with adaptive RKF45 and initial step 0.1. The exact solution is y = t3. Use an absolute tolerance of 10−8. Both embedded formulas integrate this derivative exactly apart from roundoff, so the error estimate is zero or very small. Code 1 handles zero error explicitly and limits step growth to a factor of two; observe this growth until the step is clipped at the final time.

    Connect accepted points with straight lines. Evaluate that plotted polyline on a grid of 501 times and compare it with the exact curve. Then repeat the integration with maximum step 0.05; add a maximum-step limit to Code 1, as in Exercise 3. Compare the error at accepted integration points with the error of the plotted lines. A library solver’s dense output provides a third comparison. Explain why adding linearly interpolated points to the first polyline cannot recover the missing curvature.

  6. A solution that becomes infinite in finite time

    Solve y′ = y2, y(0) = 1, whose exact solution is y = 1/(1 − t) for t < 1. Integrate separately to t = 0.5, 0.9, and 0.99.

    Use absolute tolerance 10−8 in the companion code; with a library, also try relative tolerance 10−8. Record errors, smallest accepted steps, and work. Examine the exact solution before attempting a run to t = 1.1. Explain why no finite solution continues through t = 1 and why a solver failure near that point need not indicate faulty code. Keep a step or trial limit when exploring failure behavior.

  7. A fast transient inside a slowly varying system

    Integrate the system on 0 ≤ t ≤ 5, with u(0) = 1, v(0) = 0:

    u′ = −u,
    v′ = −1000(v − u) − u.

    The exact solution is u = e−t, v = e−t − e−1000t. Compare Code 2 with a tested stiff library solver listed on the resources page. Use a common absolute tolerance of 10−8, setting the library’s relative tolerance to zero if supported, or documenting its actual setting. Check errors and work, particularly after t = 0.05. For an implicit solver, supply the constant Jacobian [(−1,0); (999,−1000)] if the interface permits. Explain why the explicit solver may remain restricted after the visible fast transient has decayed.

Physics projects

The first six projects use familiar classical models; the last two offer more advanced work in rigid-body rotation and quantum dynamics. For each project, state assumptions and units, identify the state vector, and write all first-order equations before programming. Present the requested plots and a short physical interpretation. Verify the calculation using the check supplied, then repeat with a smaller step or tighter tolerances. Report the solver settings and numerical uncertainty in the quantities you extract.

Adaptive RKF45 is a useful starting point for these nonstiff default parameters. Use fixed-step velocity Verlet for the coupled masses if desired. Variables with different physical units need appropriate scaling. Limit the maximum step or use dense output for smooth plots; locate turning points and switches more carefully than simply choosing the nearest saved point.

1. A radioactive decay chain

Objective: follow the buildup and decay of an intermediate daughter species, rather than a single exponential decay.

NA′ = −λANA,
NB′ = λANA − λBNB,
NC′ = λBNB.

Assume each decay produces one daughter and the final species C is stable. Normalize the initial number of nuclei to one: (NA,NB,NC)(0) = (1,0,0). Use λA = 0.2 s−1, λB = 0.5 s−1, and integrate to 30 s. Start from Code 2 with three equations. Plot the three fractions and the two activities λANA and λBNB. Locate the time of maximum daughter abundance using its derivative or a refined search near the sampled maximum.

Validation: Check nonnegative fractions and NA + NB + NC = 1. For unequal decay constants, NA = e−λAt and NB = λA(e−λAt − e−λBt)/(λB − λA). The daughter maximum occurs at ln(λB/λA)/(λB − λA) ≈ 3.054302 s.

Possible extension: Set the two decay constants equal. Derive and test the limit NB = λt e−λt, rather than evaluating a formula with a zero denominator.

2. Thermal exchange between two finite bodies

Objective: model the evolution of both temperatures when neither body acts as a fixed-temperature reservoir.

C1T1′ = −G(T1 − T2),
C2T2′ = G(T1 − T2).

The bodies have constant heat capacities C1 = 500 J/K and C2 = 1000 J/K, with thermal conductance G = 5 W/K. Assume uniform temperature within each body and no heat exchange with the surroundings. Start at T1 = 400 K, T2 = 300 K, and integrate to 500 s. Use Code 2 with two equations. Plot both temperatures and the heat flow. Vary the heat-capacity ratio and explain which body changes temperature more.

Validation: The conserved quantity is C1T1 + C2T2, with an arbitrary common energy reference. Define Teq = (C1T1(0) + C2T2(0))/(C1 + C2) and τ = 1/[G(1/C1 + 1/C2)]. Then D = 100 K exp(−t/τ), T1 = Teq + C2D/(C1 + C2), and T2 = Teq − C1D/(C1 + C2). The defaults give Teq = 333.333333 K and τ = 66.666667 s.

Possible extension: Replace G by the positive conductance G0[1 + ((T1 − T2)/(100 K))2]. Does the final equilibrium change, or only the route to it?

3. An RLC circuit responding to a voltage pulse

Objective: follow energy supplied, stored, and dissipated during switching.

q′ = I,
L I′ = V(t) − RI − q/C.

Use a series circuit with inductance L = 0.10 H, capacitance C = 1.0 × 10−3 F, and resistance R = 5 Ω. Apply V(t) = V0 = 10 V from t = 0 to 0.10 s, then zero voltage through t = 1 s. Initially q = I = 0. Modify Code 3, treating charge and current as the two state variables. Integrate to the switching time and restart with the same state. Plot current, capacitor voltage, and stored energy. Repeat with R = 20 and 40 Ω; compare overshoot and settling behavior. Use output spacing or a maximum step small enough to resolve the oscillations.

Validation: With E = LI2/2 + q2/(2C), verify dE/dt = VI − RI2, preferably through its integrated energy balance. Before switching, the underdamped default has q = CV0[1 − e−αt(cos(ωdt) + (α/ωd)sin(ωdt))], where α = R/(2L) and ωd = sqrt(1/(LC) − α2). Critical damping occurs at R = 2 sqrt(L/C) = 20 Ω.

Possible extension: Vary pulse duration and compare the energy left in the circuit at switch-off. Keep charge and current continuous at each switch.

4. Energy transfer between two coupled masses

Objective: observe beating and energy exchange between coupled oscillators.

x1′ = v1,   x2′ = v2,
m v1′ = −kx1 − kc(x1 − x2),
m v2′ = −kx2 − kc(x2 − x1).

Two equal masses are connected to fixed supports by equal springs and to each other by a coupling spring. Displacements are measured from equilibrium. Use m = 1 kg, k = 1 N/m, kc = 0.10 N/m, x1(0) = 0.10 m, and all other initial state variables zero. Integrate to 200 s using Code 2 with four equations, or adapt Code 5 to two position coordinates. Plot displacements and the individual quantities Ei = mvi2/2 + kxi2/2. These individual energies exclude the coupling spring; show its energy separately. Relate the slow transfer to the difference between the two normal-mode frequencies.

Validation: The exact displacements are x1 = (0.10 m/2)[cos(ω+t) + cos(ω−t)] and x2 = (0.10 m/2)[cos(ω+t) − cos(ω−t)], where ω+ = sqrt(k/m) is the in-phase mode frequency and ω− = sqrt((k + 2kc)/m) is the out-of-phase mode frequency. The labels + and − describe the relative signs of the two displacements, not the ordering of the frequencies. Verify conservation of E1 + E2 + kc(x1 − x2)2/2.

Possible extension: Make the masses unequal. How does detuning change the extent of energy transfer? Derive the modified acceleration equations before changing parameters.

5. A charged particle in crossed electric and magnetic fields

Objective: separate rapid cyclotron motion from the drift of the guiding center.

x′ = vx,   y′ = vy,
vx′ = E + vy,   vy′ = −vx.

Use reduced units with positive charge q = 1, mass m = 1, magnetic field Bz = 1, and electric field Ex = E = 0.2. All initial positions and velocities are zero. Neglect radiation and relativity. Integrate to 20π using Code 2 or Code 4. Plot the trajectory and velocity components, then average the velocity over complete cyclotron periods. Explain why a magnetic field changes direction but does no work. The force depends on velocity, so the standard position-force update in Code 5 does not apply directly.

Validation: For these initial data, vx = E sin t, vy = E(cos t − 1), x = E(1 − cos t), and y = E(sin t − t). The mean drift is (0,−E). Check conservation of (vx2 + vy2)/2 − Ex, which initially equals zero.

Possible extension: Restore arbitrary q, m, and B in the equations and reverse the charge. Compare cyclotron rotation with the drift E × B/|B|2, which is independent of charge sign.

6. A rocket that loses mass during a finite burn

Objective: model propulsion with changing mass and connect the powered and coasting stages.

h′ = v,   m′ = −μ,
v′ = uμ/m − g   (during the burn).

Assume a vertical flight with constant gravity, no air resistance, and effective exhaust speed u = 2000 m/s. Use initial mass m0 = 1000 kg, dry mass md = 600 kg, positive mass-loss rate μ = 5 kg/s, and g = 9.81 m/s2. Start with h = v = 0. The thrust uμ = 10000 N exceeds the initial weight. Integrate the three equations to burnout at tb = (m0 − md)/μ = 80 s. Then set thrust and mass loss to zero and follow the coast to the apex. Plot mass, velocity, altitude, and acceleration. Use Code 2 with three components and sensible component scales or tolerances.

Validation: During the burn, m = m0 − μt and v = u ln(m0/m) − gt. Check this velocity curve. After burnout, v = vb − g(t − tb), so the additional height is vb2/(2g). Mass, altitude, and velocity remain continuous at burnout; acceleration changes abruptly. Do not let the mass decrease below the dry mass. The default parameters give vb ≈ 236.851 m/s, hb ≈ 6009.850 m, and an apex at h ≈ 8869.102 m and t ≈ 104.144 s, measured from launch.

Possible extension: Keep the fuel mass and exhaust speed fixed, then vary the burn rate while retaining immediate lift-off. Explain the effect of burn duration on gravitational velocity loss. The thrust model is discussed in NASA’s rocket thrust notes.

7. Torque-free rotation and the intermediate-axis instability

Objective: test the stability of rotation about the three principal axes of an asymmetric rigid body.

I1ω1′ = (I2 − I3)ω2ω3,
I2ω2′ = (I3 − I1)ω3ω1,
I3ω3′ = (I1 − I2)ω1ω2.

These are Euler’s equations in the body’s principal-axis frame, with no external torque. Use reduced units with principal moments (I1,I2,I3) = (1,2,3). Run three initial angular velocities: (1,0.01,0.01), (0.01,1,0.01), and (0.01,0.01,1), each to t = 50. Use Code 2 with three equations. Plot all angular-velocity components and compare the growth of the small transverse components. Reduce the perturbations to 10−3 and repeat. The components are expressed in a rotating body frame, so they need not be constant even though angular momentum in space is conserved.

Validation: Check both E = (Σ Iiωi2)/2 and |L|2 = Σ Ii2ωi2. Exact rotation about any single principal axis stays constant. Linearization near the intermediate axis predicts an early growth rate |ω2| sqrt[(I2 − I1)(I3 − I2)/(I1I3)]. Fit growth only while the transverse components remain small.

Possible extension: Compare with an axisymmetric body, such as I1 = I2 = 1, I3 = 2. Further background is available in MIT’s notes on Euler’s equations.

8. Population transfer in a two-level quantum system

Objective: solve a time-dependent quantum state as four real first-order equations and test probability conservation. This project assumes introductory quantum mechanics.

i c1′ = −(Δ/2)c1 + (Ω/2)c2,
i c2′ = (Ω/2)c1 + (Δ/2)c2.

The amplitudes describe a closed, coherently coupled two-state model. The real constants Ω and Δ have units of angular frequency and represent coupling and detuning; Planck’s constant has been divided out. Set Ω = 1 s−1 and compare Δ = 0, 1, 3 s−1. Start with c1 = 1, c2 = 0, and integrate to 20 s. To use Code 2 without complex arithmetic, write c1 = a + ib and c2 = d + ie:

a′ = −Δb/2 + Ωe/2,
b′ = Δa/2 − Ωd/2,
d′ = Ωb/2 + Δe/2,
e′ = −Ωa/2 − Δd/2.

The initial real state is (a,b,d,e) = (1,0,0,0). Keep the imaginary unit i distinct from the real component e.

Validation: Check P1 = a2 + b2, P2 = d2 + e2, and P1 + P2 = 1. The exact transferred population is P2(t) = [Ω2/(Ω2 + Δ2)] sin2(sqrt(Ω2 + Δ2)t/2). At zero detuning the first complete transfer occurs at t = π/Ω. Plot both populations and explain how detuning changes the maximum transfer.

Possible extension: Apply a finite coupling pulse by turning Ω off at a chosen time, carrying the amplitudes continuously across the switch. Compare resonant pulse durations π/(2Ω) and π/Ω. See MIT’s two-level-system notes for background.

← Back to Chapter 9