Chapter 24 Exercises — Partial Differential Equations and Finite Differences

DataField.Dev

Chapter 24 Exercises — Partial Differential Equations and Finite Differences

These exercises take you from reading a stencil to building a stable, boundary-aware heat solver you can trust. They are the point at which the running project stops being scaffolding and starts being physics, so several of them extend the solver directly — treat those as real engineering, not drills.

The difficulty tiers are:

⭐ Foundational — one idea, short code, predict-then-run.
⭐⭐ Applied — combine the stencil, stepping, stability, and boundaries into working code.
⭐⭐⭐ Challenge — a derivation, a higher-order scheme, or an implicit method; expect to think.

Problems marked † have full worked solutions in appendices/answers-to-selected.md (and, for code, in this chapter's code/exercise-solutions.f90). Odd-numbered problems are also solved there. Everything compiles with gfortran -std=f2018 -Wall. Two standing rules from the whole book apply with special force here, because a PDE solver punishes carelessness: predict the output before you run it, and never trust a stencil you have not hand-checked once — the stencil that looks right and quietly drops a factor of $h^2$ is the one that costs a research group a month.

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

A1 ⭐† Type in example-02-ftcs-1d.f90 (the 1D rod). Before compiling, predict u after one step, and write down your prediction. Then compile and run. Now change r to 0.6_dp (above the 1D limit of $1/2$), run for 40 steps, and describe — in one sentence, with the amplification factor to back it — what the printed numbers do and why.

A2 ⭐ Predict the discrete Laplacian of a linear field u(i,j) = 2.0_dp*real(i,dp) + 3.0_dp*real(j,dp) on a $5\times5$ grid with the five-point stencil. State your prediction and the reason before writing any code, then verify. (Hint: what is $\nabla^2$ of a linear function, exactly?)

A3 ⭐⭐† The snippet below runs the 2D solver for 100 steps and prints the maximum interior temperature. Predict, roughly, what it prints for r = 0.20 versus r = 0.30, and explain both. Which of the two runs would you put in a paper, and why is the other one worthless?

real(dp) :: hot
! ... set up plate, hot top edge = 100, choose r via alpha, dt, dx ...
do step = 1, 100
  call step_field(u)               ! one FTCS step, r baked in
end do
hot = maxval(u(2:n-1, 2:n-1))
print '(a, es12.3)', 'max interior T = ', hot

A4 ⭐⭐ Take the Project Checkpoint's $5\times5$ plate, but hold the left edge at 100 (instead of the top) and all others at 0. Hand-compute the interior after one step with $r = 0.2$, then confirm with a modified project-checkpoint.f90. Which interior cells warm, and does the symmetry of your answer match the symmetry of the boundary?

Part B — The Stencil and Discretization

B5 ⭐† Write the 1D three-point Laplacian $u_{i-1} - 2u_i + u_{i+1}$ over the interior of a rank-1 array u(n) as a single array-section statement (no loop), scaled by $1/\Delta x^2$, storing it in lap(2:n-1). Test it on u = (i-1)**2 and confirm every interior value is the exact $d^2(x^2)/dx^2 = 2$.

B6 ⭐⭐ Our laplacian(u, dx, dy) already handles $\Delta x \neq \Delta y$. Build a $6\times6$ field u(i,j) = real(i,dp)**2 (it varies only along the first index), give it dx = 0.5, dy = 2.0, and predict the interior Laplacian by hand before running. Explain which of the two direction-terms survives and why the answer is $8$.

B7 ⭐⭐† (Verifies the Chapter 22 project increment.) Numerically confirm the five-point stencil is second-order accurate. Take $u(x,y) = \sin(x)\sin(y)$, whose exact Laplacian is $-2\sin(x)\sin(y)$. Compute the maximum error of the discrete Laplacian on grids with $h$, $h/2$, and $h/4$, and show the error falls by about $4\times$ each time you halve $h$. Report the three errors and the two ratios, and state which number — the error or the ratio — is the meaningful result.

B8 ⭐⭐⭐ The nine-point stencil adds the four diagonal neighbours and is fourth-order accurate for the Laplacian. Look up its weights, implement it as array sections (the diagonals are slices like u(1:n-2,1:n-2)), and repeat B7's convergence study. Does the error now fall by $\sim16\times$ per halving? What did you trade for the higher order?

Part C — Stability and the CFL Condition

C9 ⭐† Compute the largest stable timestep $\Delta t$ for $\alpha = 1.0\times10^{-4}\ \mathrm{m^2/s}$ and $h = 0.01\ \mathrm{m}$, in 1D ($r \le 1/2$), 2D ($r \le 1/4$), and 3D ($r \le 1/6$). Show the numbers, and say in one sentence what the trend across dimensions costs you.

C10 ⭐⭐† Find the bug. This code "worked fine on the coarse grid" but produces NaN on the fine one:

real(dp), parameter :: alpha = 0.1_dp, dt = 0.002_dp
real(dp) :: h
h = 1.0_dp / real(n - 1, dp)      ! n was 26 last week, now 51
! ... FTCS step with r = alpha*dt/h**2 ...

Compute $r$ for n = 26 and n = 51, diagnose the failure in terms of the CFL condition, then rewrite the two lines so the code is stable on any n.

C11 ⭐⭐ Take example-03-cfl-blowup.f90 and, without running it, hand-predict the interior values at steps 1, 2, and 3 for r = 0.25 (exactly the limit) instead of 0.5. Is the checkerboard growing, shrinking, or holding? What does that tell you about running a production job exactly at the CFL limit rather than a safe fraction below it?

C12 ⭐⭐⭐† Do the von Neumann analysis for the 1D explicit scheme. Substitute $u^n_j = G^n e^{\mathrm i k j h}$ into $u^{n+1}_j = u^n_j + r(u^n_{j+1} - 2u^n_j + u^n_{j-1})$, use $1 - \cos\theta = 2\sin^2(\theta/2)$ to show the amplification factor is $G = 1 - 4r\sin^2(kh/2)$, and derive the stability limit $r \le 1/2$ from $|G| \le 1$. Then say, in one line, how the 2D version changes the result.

C13 ⭐⭐ Back-of-envelope stability instrument. Write a one-line diagnostic that, given alpha, dt, dx, dy, prints the diffusion number and a verdict (STABLE/UNSTABLE) against the 2D limit. Where in a real program should this live, and why is printing $r$ at startup worth the two lines it costs?

Part D — Boundary Conditions

D14 ⭐† Modify the 1D solver so both ends are zero-Neumann (insulated): u(1) = u(2) and u(n) = u(n-1), re-imposed each step. Start from a non-uniform interior of your choice, predict the steady state, and explain the prediction from conservation of heat. (With no heat entering or leaving, what must the final field be?)

D15 ⭐⭐ Implement a periodic 1D heat step using modulo for the neighbour indices, so the interior update wraps end-to-end. Seed a single hot node and describe how the profile evolves differently than with absorbing Dirichlet ends. What quantity is conserved here that is not conserved under Dirichlet?

D16 ⭐⭐† Design it — mixed boundaries. A plate is insulated on its left edge (zero-Neumann) and held at fixed temperatures on the other three (Dirichlet). Write the boundary-update routine that imposes all four conditions each step, and state precisely where in the time loop it must be called relative to the interior update — and what goes wrong if you call it in the other order.

D17 ⭐⭐ Port it. Here is a NumPy zero-Neumann edge update a colleague wrote:

u[0, :] = u[1, :]        # top edge insulated
u[:, 0] = u[:, 1]        # left edge insulated

Translate it to Fortran on a field_t's u(:,:), mind the 1-based indexing and the row/column-major difference, and say which NumPy axis corresponds to which Fortran index.

Part E — Design It (extend the solver)

E18 ⭐⭐† Add a run_to_steady(field, alpha, dt, tol, max_steps, nsteps) routine that calls step until the maximum change between successive fields drops below tol (or max_steps is hit), returning the number of steps in nsteps. Use it on a small Dirichlet plate and report the step count. Why is the max_steps cap not optional in defensive code?

E19 ⭐⭐ Add a constant heat source $s$ to the interior update: $u^{n+1} = u^n + \Delta t\,(\alpha\nabla^2 u + s)$. Show exactly where s enters the step routine. Does it change the CFL limit? Explain why the stability limit depends on the stencil, not the source.

E20 ⭐⭐⭐ Implement an implicit 1D heat step (backward-Euler): the update $(I - r\,L)\,\mathbf{u}^{n+1} = \mathbf{u}^n$, where $L$ is the tridiagonal second-difference operator. Build the matrix, apply Dirichlet rows, and solve each step by calling LAPACK's dgesv (Chapter 21). Verify it stays stable at r = 5.0 — twenty times the explicit limit — and state what you paid for that freedom.

E21 ⭐⭐ Give step an optional source argument (a scalar real(dp)) that adds a uniform source as in E19 when present and defaults to none when absent. Show the signature, the one-line present guard, and a call site that uses it — reusing the optional-argument machinery from Chapter 6.

Part F — Modernize It

F22 ⭐⭐† The fragment below is the heart of a real FORTRAN 77 Jacobi relaxation kernel for steady-state heat — exactly the five-point stencil, wearing 1977 clothes. Modernize it to the book's style (implicit none, free-form, real(dp), array section, intent), and identify which of §24's boundary conditions the DO bounds are silently assuming.

      DO 20 J = 2, N-1
      DO 10 I = 2, M-1
      UNEW(I,J) = 0.25*(U(I-1,J)+U(I+1,J)+U(I,J-1)+U(I,J+1))
   10 CONTINUE
   20 CONTINUE

F23 ⭐⭐ The relaxation kernel of F22 has no explicit timestep and no $\alpha$ — the new value is just the average of the four neighbours. Explain how that is the $r = 1/4$ FTCS update in disguise (substitute $r = 1/4$ into the 2D update and simplify), and why the old-timers picked exactly $1/4$: what does relaxation converge to, and why is $1/4$ the fastest safe choice?

Part G — Back of the Envelope

G24 ⭐† You must simulate $T = 10$ seconds of physical diffusion with $\alpha = 10^{-4}\ \mathrm{m^2/s}$ on a $1\ \mathrm{m}$ plate gridded at $h = 0.005\ \mathrm{m}$. Using a CFL-safe 2D timestep, roughly how many time steps is that? (Order of magnitude is fine; show the arithmetic.)

G25 ⭐⭐ A $10{,}000 \times 10{,}000$ real(dp) field is how many gigabytes? The FTCS step needs the current field plus one work array. What is the minimum memory footprint, and what does that imply about a laptop versus a cluster node — and about Chapter 34's domain decomposition?

G26 ⭐⭐† You refine a 2D explicit heat run from spacing $h$ to $h/2$ for a sharper picture, keeping the physics and total simulated time fixed. By what factor does the total work (grid points $\times$ time steps) increase? Work the 1D and 3D cases too, express the general rule as a power of $h$, and explain in one sentence why explicit diffusion codes dread fine grids.

Part H — Interleaved (earlier chapters)

H27 ⭐† (Chapter 5.) A colleague's stencil loop is do i = 2, n-1; do j = 2, n-1; u_new(i,j) = …. On a large grid it is several times slower than it should be. Name the problem in one phrase, give the fixed loop nesting, and justify it from Fortran's memory order.

H28 ⭐⭐ (Chapter 9.) The Chapter 6 step took a bare field(:,:); this chapter's takes a type(field_t). Concretely, what does the field_t version have access to that the bare-array version did not, and why does the real stencil (unlike the Chapter 6 placeholder) actually need it?

H29 ⭐⭐† (Chapter 23.) Explain, in three or four sentences, why the explicit heat step is exactly forward-Euler applied to the method-of-lines ODE system. Then say what changes — and what stays the same — if you swap forward-Euler for RK4.

H30 ⭐⭐⭐† (Chapter 13.) Add defensive validation to a configure routine: compute $r = \alpha\Delta t/h^2$ and error stop with a helpful message if $r > 0.25$ (2D). Explain why catching this at setup is far better than discovering it as NaN a thousand steps into an overnight run, and why error stop (not stop) is the right choice for a batch job.

Solutions to the starred, odd-numbered, and †-marked problems are in appendices/answers-to-selected.md; the compilable ones are in code/exercise-solutions.f90. If your output disagrees with a hand-computed "Expected output," you have found either your bug or ours — either way, you have learned the stencil.