Chapter 23 Exercises: Ordinary Differential Equations

These problems build fluency with the three things that matter: writing a right-hand side, choosing an integrator, and reasoning about its error. Work them with a compiler open. For every "predict the output" problem, write your prediction down before you compile — the discipline of hand-tracing is the skill.

Difficulty tiers. ⭐ warm-up (minutes); ⭐⭐ core (should take real thought); ⭐⭐⭐ extension (open-ended, closer to real work).

Solutions policy. Problems marked † and all odd-numbered problems have full solutions in appendices/answers-to-selected.md; the coding ones are also provided compilable in code/exercise-solutions.f90. Everything else is left for you (or your instructor).

Reuse the euler_step, rk4_step, and rk4_sys routines from the chapter's code/ directory rather than retyping them — passing the RHS as a procedure argument is exactly what lets you.

Part A — Type, Compile, and Run (predict first)

A1 (⭐) † Take example-01-euler.f90 and change the step to h = 0.5, integrating $y' = y$, $y(0) = 1$ to $t = 1$ (two steps). Predict $y_1$ and $y_2$ by hand, then run it. By how much does the larger step worsen the final error compared with the chapter's $h = 0.25$ run?

A2 (⭐) Change the RHS in example-01-euler.f90 to $f(t, y) = -2y$ (decay), keep $y(0) = 1$ and $h = 0.25$, and take one step. Predict $y_1$. Then predict what a second step gives, and say whether Euler is over- or under-shooting the true $y(t) = e^{-2t}$.

A3 (⭐) † Using rk4_step, take a single RK4 step of $y' = y$ from $y(0) = 1$ with $h = 0.5$. Predict all four stages $k_1, k_2, k_3, k_4$ and the result $y(0.5)$ by hand, then confirm against $e^{0.5} = 1.6487213$. (Answer check: the stages are $1,\ 1.25,\ 1.3125,\ 1.65625$.)

A4 (⭐) Implement the midpoint method (RK2) and integrate $y' = y$ to $t = 1$ with $h = 1$ (one step). Predict the result, then compare the ladder Euler ($2.0$) → midpoint (?) → RK4 ($2.708$) → exact ($2.718$).

A5 (⭐⭐) Modify example-02-rk4.f90 to integrate $y' = y$ to $t = 1$ with $h = 0.5$ (two steps) and then $h = 0.25$ (four steps). Record the error each time. What is the error ratio, and does it match the $16\times$ you expect from fourth order? Explain any discrepancy at these large steps.

Part B — Port It (translate to Fortran, then compare)

B1 (⭐⭐) † Port this Python logistic-growth integrator to Fortran, using euler_step with a logistic RHS ($r = 1$, $K = 10$, $N_0 = 1$, $h = 0.5$). Predict $N_1$ and $N_2$.

r, K, N, h = 1.0, 10.0, 1.0, 0.5
for _ in range(2):
    N = N + h * (r * N * (1 - N / K))
    print(N)

B2 (⭐⭐) The two-body problem in Python below uses a hand-written RK4 on a 4-vector state $(x, y, v_x, v_y)$. Port the RHS to a Fortran rhs_sys function and drive it with rk4_sys. You need not reproduce the plot — just get one step to agree.

def f(t, s):           # s = [x, y, vx, vy], GM = 1
    x, y, vx, vy = s
    r3 = (x*x + y*y)**1.5
    return np.array([vx, vy, -x/r3, -y/r3])

B3 (⭐⭐) MATLAB's [t, y] = ode45(@f, [0 10], y0) hides an adaptive Dormand-Prince integrator behind one line. Sketch (in prose and a signature) how you would structure the equivalent in Fortran: what does @f become, what does the [0 10] span become, and where does the tolerance live? You are describing an interface, not writing the solver.

Part C — Find the Bug

C1 (⭐⭐) † This "RK4" gives second-order accuracy, not fourth. Find the one-character bug and explain how you would have caught it without being told.

k1 = f(t,            y)
k2 = f(t + 0.5_dp*h, y + 0.5_dp*h*k1)
k3 = f(t + 0.5_dp*h, y + 0.5_dp*h*k1)   ! <-- ?
k4 = f(t + h,        y + h*k3)
yn = y + (h/6.0_dp) * (k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)

C2 (⭐⭐) This Euler loop updates time in the wrong place. Predict the wrong trajectory it produces for a non-autonomous RHS like $f(t, y) = t - y$, and fix it.

do n = 1, nsteps
  t = t + h
  y = y + h * f(t, y)     ! <-- t has already been advanced
end do

C3 (⭐⭐) † A student writes a system RHS whose result is declared real(dp) :: dydt(2) (fixed size) and calls rk4_sys with a 4-component orbital state. What happens — a compile error, a run-time crash, or silently wrong numbers? Explain, and give the assumed-shape fix.

C4 (⭐⭐) This code integrates a stiff decay $y' = -1000\,y$ with RK4 and h = 0.01, and the output is NaN. The coefficients are correct. Diagnose the real problem (hint: it is not a typo) and state two fixes.

Part D — Modernize It

D1 (⭐⭐) † Modernize this FORTRAN 77 Euler integrator: give it implicit none, free-form layout, a module, intent on every argument, and a dp kind. Pass the RHS as a procedure argument instead of hard-coding it.

      SUBROUTINE EULER(Y0, H, N, YOUT)
      REAL Y0, H, YOUT, Y
      INTEGER N, I
      Y = Y0
      DO 10 I = 1, N
        Y = Y + H * (Y)
   10 CONTINUE
      YOUT = Y
      RETURN
      END

D2 (⭐⭐) The F77 fragment below uses a statement function for the RHS. Convert it to a modern internal or module procedure and explain why the modern form is safer.

      F(T, Y) = -0.5 * Y + T
      Y = Y + H * F(T, Y)

Part E — Design It (extend the heat solver)

E1 (⭐⭐⭐) † Replace the explicit Euler time-step in project-checkpoint.f90 with an RK4 step (rk4_sys) applied to the same heat_rhs, and take one step of $h = 0.1$ from a cold interior. Predict the three interior values and compare with the Euler step $(0.1, 0, 0)$. Why does RK4 already spread heat to points 2 and 3 in a single step when Euler does not?

E2 (⭐⭐⭐) Generalize the method-of-lines RHS to $N$ interior points on a rod (an assumed-shape array and a loop or a whole-array second difference), with configurable hot/cold boundaries. Verify it against the 3-point hand trace for $N = 3$, then run $N = 20$ and describe the profile after 50 Euler steps.

E3 (⭐⭐⭐) Instrument the MOL heat stepper to detect instability: after each step, if any $|u_i|$ exceeds a threshold (say $10^{3}$), stop with error stop and report the step number. Then set $\Delta t = 0.6 > \Delta x^2/(2\alpha) = 0.5$ and confirm your guard fires. This is the CFL limit of Chapter 24, discovered experimentally.

Part F — Back of the Envelope

F1 (⭐⭐) † You need the solution of a smooth IVP to an accuracy of $10^{-6}$ over $t \in [0, 10]$. Estimate how many steps Euler needs versus RK4, using error $\sim h$ and error $\sim h^4$ respectively (assume the error constant is $\sim 1$). Then estimate the ratio of total RHS evaluations, remembering RK4 costs four per step.

F2 (⭐⭐) A stiff mode has $\lambda = -10^{6}$. How small must explicit Euler's step be for stability? How many steps does that force over $t \in [0, 1]$? Contrast with backward Euler, which is stable at any step — if accuracy allows $h = 10^{-2}$, how many steps does it take?

F3 (⭐⭐) † You RK4-integrate a method-of-lines system with $m = 10^{6}$ grid points. How many RHS evaluations and how many full-length array temporaries ($k_1$–$k_4$ plus the stage arguments) does one step touch? Estimate the memory traffic per step in bytes (double precision). Is this compute- or memory-bound? (Foreshadow: Chapter 28 measures exactly this.)

F4 (⭐⭐) Halving $\Delta x$ in an explicit heat solver forces $\Delta t$ to shrink by the CFL rule $\Delta t \le \Delta x^2/(2\alpha)$. If you double the spatial resolution, by what factor does the total work to reach a fixed final time increase? (Count both the extra points and the extra steps.)

Part G — Interleaved (mix earlier chapters)

G1 (⭐⭐) † (Chapter 20) You integrate for $10^{6}$ steps. Explain why the state must be real(dp) and not default real: estimate the accumulated round-off in each and argue which is tolerable. (Recall Chapter 20's error-budget reasoning.)

G2 (⭐⭐) (Chapter 6) Rewrite the scalar grow RHS as a pure function and explain, using Chapter 6's rules, exactly what pure promises and why an ODE RHS is a natural candidate for it. What would break the pure guarantee?

G3 (⭐⭐⭐) † (Chapter 21) Write out the backward-Euler step for the 3-point heat MOL as a linear system $(I - \Delta t\,A)\,\mathbf{u}^{n+1} = \mathbf{u}^{n}$. Form the $3\times3$ matrix explicitly for $\alpha = \Delta x = 1$, $\Delta t = 0.1$ (include the boundary contributions in the right-hand side), and describe how you would solve it each step with LAPACK's dgesv (or the tridiagonal dgtsv).

G4 (⭐⭐) (Chapter 22) Show that RK4's weights $\tfrac16, \tfrac13, \tfrac13, \tfrac16$ are exactly Simpson's rule applied to $\int_{t_n}^{t_n+h} f\,dt$. Why should the weights of a time integrator be quadrature weights at all?

G5 (⭐⭐⭐) (Chapter 5) Rewrite the $N$-point method-of-lines RHS of E2 with no explicit loop — a single whole-array expression using array sections u(1:n-2), u(2:n-1), u(3:n). Confirm it gives the same numbers as the loop version, and say why the array form is both clearer and friendlier to the optimizer.

Full solutions to the odd-numbered and † problems: appendices/answers-to-selected.md; the coding ones are compilable in code/exercise-solutions.f90.