HW 3 Solutions

Published

October 1, 2026

Overview

Load Environment

The following code loads the environment and makes sure all needed packages are installed. This should be at the start of most Julia scripts.

import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using Plots
using LaTeXStrings

Problems (Total: 50 Points)

Problem 1 (15)

Problem 1.1 (4)

\[\frac{dC}{dt} = 5 - 0.25C\]

(a) Denoting the right-hand side \(f(C) = 5 - 0.25C\), the forward Euler update is \[C(t + \Delta t) = C(t) + \Delta t\left(5 - 0.25\,C(t)\right).\]

(b) This is a linear system with a single rate constant, so its characteristic timescale is \[1/|-0.25| = 4\]. Since \(\Delta t = 1\) is four times smaller than the timescale, it is a reasonable starting choice, though it might not be small enough to account for discretization error.

(c) Assuming first-order convergence (as expected from the Taylor approximation underlying Forward Euler discretization), error scales linearly with \(\Delta t\). So halving \(\Delta t\) from 2 to 1 should halve the error: \[\text{error}(\Delta t = 1) \approx 0.6 \times \frac{1}{2} = 0.3.\]

Problem 1.2 (5)

\[\frac{dC}{dt} = \frac{Q}{V}\left(C_\text{in} - C\right) - kC\]

(a) With \(Q/V = 50/200 = 0.25\ \text{yr}^{-1}\), the forward Euler update is \[\begin{aligned} C(t + \Delta t) &= C(t) + \Delta t\left[\frac{Q}{V}\left(C_\text{in} - C(t)\right) - kC(t)\right] \\ &= C(t) + \Delta t\left[0.25\left(2 - C(t)\right) - 0.3\,C(t)\right] \\ &= C(t) + \Delta t\left(0.5 - 0.55 C(t)\right). \end{aligned}\]

(b) \(\frac{dC}{dt} = 0.5 - 0.55C\): the two loss terms combine into a single effective rate \(-0.55\ \text{yr}^{-1}\). The characteristic timescale is \(1/|-0.55| \approx 1.82\ \text{yr}\). \(\Delta t = 1\ \text{yr}\) could be a reasonable starting choice, though as with Problem 1.1, there could still be discretization error that would need to be reduced by refining \(\Delta t\).

(c) Here \(\Delta t\) shrinks by a factor of 4 (from 1.6 to 0.4), so under first-order convergence the error should also shrink by a factor of 4: \[\text{error}(\Delta t = 0.4) \approx 0.8 \times \frac{1}{4} = 0.2.\]

Problem 1.3 (6)

\[\frac{dx}{dt} = 3 - 0.5x^2, \qquad x \geq 0\]

(a) Denoting \(f(x) = 3 - 0.5x^2\), the forward Euler update is \[x(t + \Delta t) = x(t) + \Delta t\left(3 - 0.5\,x(t)^2\right).\]

(b) The equilibrium solves \(3 - 0.5(x^*)^2 = 0 \implies x^* = \sqrt{6} \approx 2.449\) (taking the positive root, per the problem’s domain). \(f'(x) = -x\), so \(g = f'(x^*) = -\sqrt{6} \approx -2.449\). The local characteristic timescale there is \(1/|\sqrt{6}| \approx 0.408\).

\(\Delta t = 0.5\) is larger than this local timescale so this is not a reasonable starting value near this equilibrium.

(c) Here \(\Delta t\) shrinks by a factor of 3, so \[\text{error}(\Delta t = 0.1) \approx 1.5 \times \frac{1}{3} = 0.5.\]

Problem 2 (35)

Problem 2.1 (5)

Denoting the right-hand sides \(f(P,M) = L - s_\text{out}P - hP + \kappa M \frac{P^q}{m^q+P^q}\) and \(g(P,M) = hP - bM - \kappa M \frac{P^q}{m^q+P^q}\), the coupled forward Euler update is

\[\begin{aligned} P(t + \Delta t) &= P(t) + \Delta t\left[L - s_\text{out}P(t) - hP(t) + \kappa M(t)\,\frac{P(t)^q}{m^q + P(t)^q}\right] \\ M(t + \Delta t) &= M(t) + \Delta t\left[hP(t) - bM(t) - \kappa M(t)\,\frac{P(t)^q}{m^q + P(t)^q}\right] \end{aligned}\]

Note: both updates must use the recycling term evaluated at the current \(P(t)\) and \(M(t)\). A common mistake is to update \(P\) first and then use the new \(P(t+\Delta t)\) when updating \(M\) — this is not forward Euler (it’s closer to a Gauss-Seidel-style update), which requires joint linearization (recall the multivariate Taylor expansion we used to derive the condition for higher-dimensional stability). As a result, it does not have the same error properties we derived in class. You won’t lose points for this — we didn’t talk about it in class, but it’s good to be aware of.

Problem 2.2 (13)

Choosing \(T\). The slow variable is the sediment pool \(M\), whose timescale is \(1/b = 1/0.02 = 50\ \text{yr}\). \(T = 300\ \text{yr}\) is six of those — comfortably long enough for \(M\) to fill.

Choosing the quantity of interest. \(P(30)\), while the lake is still changing. The value at the end of the run would be useless here: by \(t = 300\) the system has relaxed, and every stable step size returns the same number to six or more digits.

outflow_rate, sediment_transfer = 0.3, 0.3
burial_rate, recycling_coefficient = 0.02, 0.15
half_saturation, steepness = 25.0, 8.0

function recycling_flux(phosphorus, sediment)
    saturation = phosphorus^steepness /
                 (half_saturation^steepness + phosphorus^steepness)
    recycling_coefficient * sediment * saturation
end

function simulate_lake(loading, horizon, step)
    steps = Int(round(horizon / step))
    trajectory = zeros(steps)
    phosphorus, sediment = 0.0, 0.0
    for i in 1:steps
        recycled = recycling_flux(phosphorus, sediment)
        change_in_phosphorus = loading - outflow_rate*phosphorus -
                               sediment_transfer*phosphorus + recycled
        change_in_sediment = sediment_transfer*phosphorus -
                             burial_rate*sediment - recycled
        phosphorus += step * change_in_phosphorus
        sediment += step * change_in_sediment
        trajectory[i] = phosphorus
    end
    trajectory
end

phosphorus_at_30(step) = simulate_lake(10.0, 300.0, step)[Int(30 / step)]

step_reference = 0.001
reference_value = phosphorus_at_30(step_reference)

step_sizes = [2.0, 1.0, 0.5, 0.25, 0.125, 0.0625]
values = [phosphorus_at_30(step) for step in step_sizes]
errors = [abs(value - reference_value) for value in values]
orders = [log2(errors[i] / errors[i+1]) for i in 1:length(errors)-1]

Reference solution at \(\Delta t = 0.001\) yr gives \(P(30) =\) 18.58479 \(\mu\text{g}/\text{L}\).

\(\Delta t\) [yr] \(P(30)\) error empirical order
2.0000 18.6349 5.01e-02 —
1.0000 18.6092 2.44e-02 1.04
0.5000 18.5970 1.22e-02 1.00
0.2500 18.5909 6.09e-03 1.00
0.1250 18.5878 3.04e-03 1.00
0.0625 18.5863 1.51e-03 1.01
Figure 1: Error in P(30) against step size. The slope is the order of convergence.

The error is first order, as expected for forward Euler — and the errors span a factor of 33 from coarsest to finest. \(\Delta t = 0.125\ \text{yr}\) is used for the rest of this problem: the error there is 0.003 \(\mu\text{g}/\text{L}\), far below anything that would change the interpretation in Problem 2.3, and the run is still fast.

NoteThe horizon and the step size are set by different processes

\(T\) was chosen based on the slowest process in the system, the sediment pool at \(1/b = 50\ \text{yr}\). \(\Delta t\) is constrained by the fastest, which is \(P\) itself at \(1/(s_\text{out} + h) \approx 1.7\ \text{yr}\). A step chosen by looking only at \(M\) would be wildly too large.

Steps well above the fast timescale do not merely lose accuracy — they can destabilise. At \(\Delta t = 4\ \text{yr}\) this model blows up entirely. Note the transition is not a perfectly abrupt: \(\Delta t = 3.5\) still returns a finite (but badly wrong) answer.

Problem 2.3 (10)

Simulating at \(L = 10\ \mu\text{g}/(\text{L}\cdot\text{yr})\) with the step size chosen in Problem 2.2 gives the trajectory below.

Figure 2: P(t) for the two-stock lake model at L = 10 μg/(L·yr).

From Problem 2.2, \(P(30) =\) 18.6 \(\mu\text{g}/\text{L}\) — a modest, plausibly-stable oligotrophic reading, a little above the naive one-stock estimate but not alarmingly so. Once the simulation has fully run, \(P\) settles to a much higher equilibrium of about 29.0 \(\mu\text{g}/\text{L}\) — solidly eutrophic, well above \(m = 25\).

Physically: at \(t = 30\), the sediment pool \(M\) hasn’t yet accumulated enough phosphorus for recycling to matter much, so the lake behaves close to a simple outflow-dominated system. As \(M\) slowly builds (over the sediment timescale, \(1/b = 50\ \text{yr}\)), the recycling term \(\kappa M \cdot P^q/(m^q+P^q)\) becomes large enough to meaningfully add phosphorus back to the water column. This triggers a rapid transition — visible as a sharp rise and an overshoot, peaking around \(P =\) 38.2 \(\mu\text{g}/\text{L}\) near \(t =\) 45 yr — before the system settles into its new, eutrophic equilibrium by roughly \(t = 80\)–\(90\ \text{yr}\).

Problem 2.4 (7)

The naive one-stock model, \(\frac{dP}{dt} = L - s_\text{out}P - hP\), has equilibrium \[P^*_\text{naive} = \frac{L}{s_\text{out}+h} = \frac{10}{0.6} \approx 16.7\ \mu\text{g}/\text{L}.\]

Since this is well below \(m = 25\ \mu\text{g}/\text{L}\), a manager using this simplified model (or, equivalently, only running a short simulation of the full model — recall \(P \approx 18.6\) at \(t=30\) in Problem 2.3) would conclude this loading is safely oligotrophic.

But the true long-run behavior of the full two-stock model shows the lake actually settles into a eutrophic state (\(P \approx 29\ \mu\text{g}/\text{L}\)) at this same loading, once the sediment pool has had time to fill. The naive model isn’t just quantitatively a bit off — it gets the qualitative outcome wrong.

References

List any external references consulted, including classmates.