{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# HW 3 Solutions\n",
        "\n",
        "## Overview\n",
        "\n",
        "### Load Environment\n",
        "\n",
        "The following code loads the environment and makes sure all needed\n",
        "packages are installed. This should be at the start of most Julia\n",
        "scripts."
      ],
      "id": "676e45f3-3f83-4491-bc4f-432fc6c16d40"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import Pkg\n",
        "Pkg.activate(@__DIR__)\n",
        "Pkg.instantiate()"
      ],
      "id": "2"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "using Plots\n",
        "using LaTeXStrings"
      ],
      "id": "4"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problems (Total: 50 Points)\n",
        "\n",
        "### Problem 1 (15)\n",
        "\n",
        "#### Problem 1.1 (4)\n",
        "\n",
        "$$\\frac{dC}{dt} = 5 - 0.25C$$\n",
        "\n",
        "**(a)** Denoting the right-hand side $f(C) = 5 - 0.25C$, the forward\n",
        "Euler update is\n",
        "$$C(t + \\Delta t) = C(t) + \\Delta t\\left(5 - 0.25\\,C(t)\\right).$$\n",
        "\n",
        "**(b)** This is a linear system with a single rate constant, so its\n",
        "characteristic timescale is $$1/|-0.25| = 4$$. Since $\\Delta t = 1$ is\n",
        "four times smaller than the timescale, it is a reasonable starting\n",
        "choice, though it might not be small enough to account for\n",
        "discretization error.\n",
        "\n",
        "**(c)** Assuming first-order convergence (as expected from the Taylor\n",
        "approximation underlying Forward Euler discretization), error scales\n",
        "linearly with $\\Delta t$. So halving $\\Delta t$ from 2 to 1 should halve\n",
        "the error:\n",
        "$$\\text{error}(\\Delta t = 1) \\approx 0.6 \\times \\frac{1}{2} = 0.3.$$\n",
        "\n",
        "#### Problem 1.2 (5)\n",
        "\n",
        "$$\\frac{dC}{dt} = \\frac{Q}{V}\\left(C_\\text{in} - C\\right) - kC$$\n",
        "\n",
        "**(a)** With $Q/V = 50/200 = 0.25\\ \\text{yr}^{-1}$, the forward Euler\n",
        "update is $$\\begin{aligned}\n",
        "C(t + \\Delta t) &= C(t) + \\Delta t\\left[\\frac{Q}{V}\\left(C_\\text{in} - C(t)\\right) - kC(t)\\right] \\\\\n",
        "&= C(t) + \\Delta t\\left[0.25\\left(2 - C(t)\\right) - 0.3\\,C(t)\\right] \\\\\n",
        "&= C(t) + \\Delta t\\left(0.5 - 0.55 C(t)\\right).\n",
        "\\end{aligned}$$\n",
        "\n",
        "**(b)** $\\frac{dC}{dt} = 0.5 - 0.55C$: the two loss terms combine into a\n",
        "single effective rate $-0.55\\ \\text{yr}^{-1}$. The characteristic\n",
        "timescale is $1/|-0.55| \\approx 1.82\\ \\text{yr}$.\n",
        "$\\Delta t = 1\\ \\text{yr}$ could be a reasonable starting choice, though\n",
        "as with Problem 1.1, there could still be discretization error that\n",
        "would need to be reduced by refining $\\Delta t$.\n",
        "\n",
        "**(c)** Here $\\Delta t$ shrinks by a factor of 4 (from 1.6 to 0.4), so\n",
        "under first-order convergence the error should also shrink by a factor\n",
        "of 4:\n",
        "$$\\text{error}(\\Delta t = 0.4) \\approx 0.8 \\times \\frac{1}{4} = 0.2.$$\n",
        "\n",
        "#### Problem 1.3 (6)\n",
        "\n",
        "$$\\frac{dx}{dt} = 3 - 0.5x^2, \\qquad x \\geq 0$$\n",
        "\n",
        "**(a)** Denoting $f(x) = 3 - 0.5x^2$, the forward Euler update is\n",
        "$$x(t + \\Delta t) = x(t) + \\Delta t\\left(3 - 0.5\\,x(t)^2\\right).$$\n",
        "\n",
        "**(b)** The equilibrium solves\n",
        "$3 - 0.5(x^*)^2 = 0 \\implies x^* = \\sqrt{6} \\approx 2.449$ (taking the\n",
        "positive root, per the problem’s domain). $f'(x) = -x$, so\n",
        "$g = f'(x^*) = -\\sqrt{6} \\approx -2.449$. The local characteristic\n",
        "timescale there is $1/|\\sqrt{6}| \\approx 0.408$.\n",
        "\n",
        "$\\Delta t = 0.5$ is *larger* than this local timescale so this is not a\n",
        "reasonable starting value near this equilibrium.\n",
        "\n",
        "**(c)** Here $\\Delta t$ shrinks by a factor of 3, so\n",
        "$$\\text{error}(\\Delta t = 0.1) \\approx 1.5 \\times \\frac{1}{3} = 0.5.$$\n",
        "\n",
        "### Problem 2 (35)\n",
        "\n",
        "#### Problem 2.1 (5)\n",
        "\n",
        "Denoting the right-hand sides\n",
        "$f(P,M) = L - s_\\text{out}P - hP + \\kappa M \\frac{P^q}{m^q+P^q}$ and\n",
        "$g(P,M) = hP - bM - \\kappa M \\frac{P^q}{m^q+P^q}$, the coupled forward\n",
        "Euler update is\n",
        "\n",
        "$$\\begin{aligned}\n",
        "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] \\\\\n",
        "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]\n",
        "\\end{aligned}$$\n",
        "\n",
        "> **Note**\n",
        ">\n",
        "> Note: **both** updates must use the recycling term evaluated at the\n",
        "> *current* $P(t)$ and $M(t)$. A common mistake is to update $P$ first\n",
        "> and then use the *new* $P(t+\\Delta t)$ when updating $M$ — this is not\n",
        "> forward Euler (it’s closer to a Gauss-Seidel-style update), which\n",
        "> requires joint linearization (recall the multivariate Taylor expansion\n",
        "> we used to derive the condition for higher-dimensional stability). As\n",
        "> a result, it does not have the same error properties we derived in\n",
        "> class. You won’t lose points for this — we didn’t talk about it in\n",
        "> class, but it’s good to be aware of.\n",
        "\n",
        "#### Problem 2.2 (13)\n",
        "\n",
        "**Choosing $T$.** The slow variable is the sediment pool $M$, whose\n",
        "timescale is $1/b = 1/0.02 = 50\\ \\text{yr}$. $T = 300\\ \\text{yr}$ is six\n",
        "of those — comfortably long enough for $M$ to fill.\n",
        "\n",
        "**Choosing the quantity of interest.** $P(30)$, while the lake is still\n",
        "changing. The value at the *end* of the run would be useless here: by\n",
        "$t = 300$ the system has relaxed, and every stable step size returns the\n",
        "same number to six or more digits."
      ],
      "id": "64c0cda0-8335-4745-b42d-021714ccaf04"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "outflow_rate, sediment_transfer = 0.3, 0.3\n",
        "burial_rate, recycling_coefficient = 0.02, 0.15\n",
        "half_saturation, steepness = 25.0, 8.0\n",
        "\n",
        "function recycling_flux(phosphorus, sediment)\n",
        "    saturation = phosphorus^steepness /\n",
        "                 (half_saturation^steepness + phosphorus^steepness)\n",
        "    recycling_coefficient * sediment * saturation\n",
        "end\n",
        "\n",
        "function simulate_lake(loading, horizon, step)\n",
        "    steps = Int(round(horizon / step))\n",
        "    trajectory = zeros(steps)\n",
        "    phosphorus, sediment = 0.0, 0.0\n",
        "    for i in 1:steps\n",
        "        recycled = recycling_flux(phosphorus, sediment)\n",
        "        change_in_phosphorus = loading - outflow_rate*phosphorus -\n",
        "                               sediment_transfer*phosphorus + recycled\n",
        "        change_in_sediment = sediment_transfer*phosphorus -\n",
        "                             burial_rate*sediment - recycled\n",
        "        phosphorus += step * change_in_phosphorus\n",
        "        sediment += step * change_in_sediment\n",
        "        trajectory[i] = phosphorus\n",
        "    end\n",
        "    trajectory\n",
        "end\n",
        "\n",
        "phosphorus_at_30(step) = simulate_lake(10.0, 300.0, step)[Int(30 / step)]\n",
        "\n",
        "step_reference = 0.001\n",
        "reference_value = phosphorus_at_30(step_reference)\n",
        "\n",
        "step_sizes = [2.0, 1.0, 0.5, 0.25, 0.125, 0.0625]\n",
        "values = [phosphorus_at_30(step) for step in step_sizes]\n",
        "errors = [abs(value - reference_value) for value in values]\n",
        "orders = [log2(errors[i] / errors[i+1]) for i in 1:length(errors)-1]"
      ],
      "id": "6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reference solution at $\\Delta t = 0.001$ yr gives $P(30) =$ 18.58479\n",
        "$\\mu\\text{g}/\\text{L}$.\n",
        "\n",
        "| $\\Delta t$ \\[yr\\] | $P(30)$ |    error | empirical order |\n",
        "|------------------:|--------:|---------:|----------------:|\n",
        "|            2.0000 | 18.6349 | 5.01e-02 |               — |\n",
        "|            1.0000 | 18.6092 | 2.44e-02 |            1.04 |\n",
        "|            0.5000 | 18.5970 | 1.22e-02 |            1.00 |\n",
        "|            0.2500 | 18.5909 | 6.09e-03 |            1.00 |\n",
        "|            0.1250 | 18.5878 | 3.04e-03 |            1.00 |\n",
        "|            0.0625 | 18.5863 | 1.51e-03 |            1.01 |"
      ],
      "id": "0d4db826-6874-4355-95d3-c6cf264c499f"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [],
      "id": "cell-fig-convergence"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The error is first order, as expected for forward Euler — and the errors\n",
        "span a factor of 33 from coarsest to finest.\n",
        "**$\\Delta t = 0.125\\ \\text{yr}$** is used for the rest of this problem:\n",
        "the error there is 0.003 $\\mu\\text{g}/\\text{L}$, far below anything that\n",
        "would change the interpretation in Problem 2.3, and the run is still\n",
        "fast.\n",
        "\n",
        "> **The horizon and the step size are set by different processes**\n",
        ">\n",
        "> $T$ was chosen based on the *slowest* process in the system, the\n",
        "> sediment pool at $1/b = 50\\ \\text{yr}$. $\\Delta t$ is constrained by\n",
        "> the *fastest*, which is $P$ itself at\n",
        "> $1/(s_\\text{out} + h) \\approx 1.7\\ \\text{yr}$. A step chosen by\n",
        "> looking only at $M$ would be wildly too large.\n",
        ">\n",
        "> Steps well above the fast timescale do not merely lose accuracy — they\n",
        "> can destabilise. At $\\Delta t = 4\\ \\text{yr}$ this model blows up\n",
        "> entirely. Note the transition is not a perfectly abrupt:\n",
        "> $\\Delta t = 3.5$ still returns a finite (but badly wrong) answer.\n",
        "\n",
        "#### Problem 2.3 (10)\n",
        "\n",
        "Simulating at $L = 10\\ \\mu\\text{g}/(\\text{L}\\cdot\\text{yr})$ with the\n",
        "step size chosen in Problem 2.2 gives the trajectory below."
      ],
      "id": "71ae5221-77d8-434a-bf89-980c115c1b3a"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [],
      "id": "cell-fig-lake-flip"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From Problem 2.2, $P(30) =$ 18.6 $\\mu\\text{g}/\\text{L}$ — a modest,\n",
        "plausibly-stable oligotrophic reading, a little above the naive\n",
        "one-stock estimate but not alarmingly so. Once the simulation has fully\n",
        "run, $P$ settles to a **much higher** equilibrium of about 29.0\n",
        "$\\mu\\text{g}/\\text{L}$ — solidly eutrophic, well above $m = 25$.\n",
        "\n",
        "Physically: at $t = 30$, the sediment pool $M$ hasn’t yet accumulated\n",
        "enough phosphorus for recycling to matter much, so the lake behaves\n",
        "close to a simple outflow-dominated system. As $M$ slowly builds (over\n",
        "the sediment timescale, $1/b = 50\\ \\text{yr}$), the recycling term\n",
        "$\\kappa M \\cdot P^q/(m^q+P^q)$ becomes large enough to meaningfully add\n",
        "phosphorus back to the water column. This triggers a rapid transition —\n",
        "visible as a sharp rise and an overshoot, peaking around $P =$ 38.2\n",
        "$\\mu\\text{g}/\\text{L}$ near $t =$ 45 yr — before the system settles into\n",
        "its new, eutrophic equilibrium by roughly $t = 80$–$90\\ \\text{yr}$.\n",
        "\n",
        "#### Problem 2.4 (7)\n",
        "\n",
        "The naive one-stock model, $\\frac{dP}{dt} = L - s_\\text{out}P - hP$, has\n",
        "equilibrium\n",
        "$$P^*_\\text{naive} = \\frac{L}{s_\\text{out}+h} = \\frac{10}{0.6} \\approx 16.7\\ \\mu\\text{g}/\\text{L}.$$\n",
        "\n",
        "Since this is well below $m = 25\\ \\mu\\text{g}/\\text{L}$, a manager using\n",
        "this simplified model (or, equivalently, only running a short simulation\n",
        "of the full model — recall $P \\approx 18.6$ at $t=30$ in Problem 2.3)\n",
        "would conclude this loading is **safely oligotrophic**.\n",
        "\n",
        "But the true long-run behavior of the full two-stock model shows the\n",
        "lake actually settles into a **eutrophic** state\n",
        "($P \\approx 29\\ \\mu\\text{g}/\\text{L}$) at this same loading, once the\n",
        "sediment pool has had time to fill. The naive model isn’t just\n",
        "quantitatively a bit off — it gets the **qualitative outcome wrong**."
      ],
      "id": "4f0db0f6-196d-436b-b8df-47a61e96ab10"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n",
        "\n",
        "List any external references consulted, including classmates."
      ],
      "id": "8a9e5a23-1874-44bb-99c3-c571f9ddf4bb"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernel_info": {
      "name": "julia"
    },
    "kernelspec": {
      "name": "julia",
      "display_name": "Julia",
      "language": "julia"
    },
    "language_info": {
      "name": "julia",
      "codemirror_mode": "julia",
      "version": "1.11.5"
    }
  }
}