{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# BEE 4750 Homework 3: Discretization, Numerical Convergence, and\n", "\n", "Simulation\n", "\n", "**Name**:\n", "\n", "**ID**:\n", "\n", "> **Due Date**\n", ">\n", "> Thursday, 09/24/26, 9:00pm\n", "\n", "## Overview\n", "\n", "### Instructions\n", "\n", "- Problem 1 asks you to discretize several systems models by hand and\n", " reason about characteristic timescales and convergence, without\n", " writing any code.\n", "- Problem 2 asks you to extend the shallow lake model from class to\n", " include a sediment phosphorus stock, then use numerical simulation to\n", " check convergence and study the resulting long-term lake dynamics.\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": "d4aa58f8-136c-4746-9476-66798a26db26" }, { "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)" ], "id": "8d6aae84-cd08-4800-ac9f-90fc8eef3c4e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1 (15)\n", "\n", "For each of the following systems models, (i) write the forward Euler\n", "update rule, (ii) find the (local) characteristic timescale and use it\n", "to evaluate whether the given step size $\\Delta t$ is a reasonable\n", "choice, and (iii) use the given error to predict the error at the new\n", "step size, assuming the method is well into its first-order convergence\n", "regime." ], "id": "117bc143-fd3f-40d2-84c5-6834bf59b102" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 1.1 (4)\n", "\n", "$$\\frac{dC}{dt} = 5 - 0.25C$$\n", "\n", "1. Write the forward Euler update rule for $C(t + \\Delta t)$.\n", "\n", "2. Find the characteristic timescale of this system. Is $\\Delta t = 1$\n", " a reasonable starting choice of step size? Justify your answer.\n", "\n", "3. A forward Euler simulation of this system has an error of $0.6$ at\n", " $\\Delta t = 2$. What error would you expect at $\\Delta t = 1$?" ], "id": "78c99635-ce64-4a82-9d71-29666c9cbc0c" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 1.2 (5)\n", "\n", "A well-mixed holding tank has volume $V = 200\\ \\text{m}^3$ and a\n", "through-flow of $Q = 50\\ \\text{m}^3/\\text{yr}$ carrying an inflow\n", "concentration $C_\\text{in} = 2\\ \\text{mg/L}$. A pollutant in the tank\n", "also decays at a first-order rate $k = 0.3\\ \\text{yr}^{-1}$. The\n", "concentration in the tank, $C(t)$, is governed by\n", "\n", "$$\\frac{dC}{dt} = \\frac{Q}{V}\\left(C_\\text{in} - C\\right) - kC.$$\n", "\n", "1. Write the forward Euler update rule for $C(t + \\Delta t)$.\n", "\n", "2. Find the characteristic timescale of this system. Is\n", " $\\Delta t = 1\\ \\text{yr}$ a reasonable starting choice of step size?\n", " Justify your answer.\n", "\n", "3. A forward Euler simulation of this system has an error of $0.8$ at\n", " $\\Delta t = 1.6\\ \\text{yr}$. What error would you expect at\n", " $\\Delta t = 0.4\\ \\text{yr}$?" ], "id": "00ac3ed5-172a-423e-bc12-afdc5d38ec12" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 1.3 (6)\n", "\n", "$$\\frac{dx}{dt} = 3 - 0.5x^2, \\qquad x \\geq 0$$\n", "\n", "1. Write the forward Euler update rule for $x(t + \\Delta t)$.\n", "\n", "2. This system’s rate of change depends on the state $x$, so there is\n", " no single characteristic timescale for the whole system — only a\n", " *local* one. Find the system’s equilibrium (be careful about the\n", " problem domain), then linearize near it to find the local feedback\n", " gain and the local characteristic timescale there. Is\n", " $\\Delta t = 0.5$ a reasonable starting choice of step size near this\n", " equilibrium? Justify your answer.\n", "\n", "3. A forward Euler simulation of this system has an error of $1.5$ at\n", " $\\Delta t = 0.3$. What error would you expect at $\\Delta t = 0.1$?" ], "id": "90e189d4-58cb-435d-a1dd-b648b4d412bc" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 2 (35)\n", "\n", "In class, we modeled a shallow lake’s phosphorus (P) concentration with\n", "\n", "$$\\frac{dP}{dt} = L - sP + R(P),$$\n", "\n", "treating the maximum sediment release rate $r$ (inside $R(P)$) as a\n", "fixed constant. A more advanced version of this model could be built by\n", "noting that $r$ is really $\\kappa M$, where $M$ is a **sediment\n", "phosphorus pool** — a slowly-filling stock.\n", "\n", "The two-stock model splits the single loss rate $s$ into two pieces\n", "(outflow/permanent burial, $s_\\text{out}$, and transfer into the\n", "sediment, $h$), and adds a second state variable $M(t)$ for the sediment\n", "pool:\n", "\n", "$$\\begin{aligned}\n", "\\frac{dP}{dt} &= L - s_\\text{out}P - hP + \\kappa M\\,\\frac{P^q}{m^q + P^q} \\\\[0.3em]\n", "\\frac{dM}{dt} &= hP - bM - \\kappa M\\,\\frac{P^q}{m^q + P^q}\n", "\\end{aligned}$$\n", "\n", "where $b$ is the rate at which sediment P is permanently buried. Use the\n", "parameter values in\n", "Table 1 for the rest\n", "of this problem.\n", "\n", "| Parameter | Meaning | Value |\n", "|:------------------------:|:--------------------|:------------------------:|\n", "| $s_\\text{out}$ | P outflow/permanent burial rate (from the water column) | $0.3\\ \\text{yr}^{-1}$ |\n", "| $h$ | P transfer rate into the sediment | $0.3\\ \\text{yr}^{-1}$ |\n", "| $b$ | sediment P burial rate | $0.02\\ \\text{yr}^{-1}$ |\n", "| $\\kappa$ | sediment recycling rate coefficient | $0.15\\ \\text{yr}^{-1}$ |\n", "| $m$ | P concentration at half-maximum sediment release | $25\\ \\mu\\text{g}/\\text{L}$ |\n", "| $q$ | steepness of the sediment release response | $8$ |\n", "\n", "Table 1: Parameters for the two-stock lake model.\n", "\n", "Assume the lake starts with no phosphorus anywhere in the system,\n", "$P(0) = M(0) = 0$." ], "id": "9827138c-ba24-4c06-9b21-a260e68b3aca" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 2.1 (5)\n", "\n", "Derive the forward Euler update rules for $P(t+\\Delta t)$ and\n", "$M(t + \\Delta t)$. Note that these two updates are coupled: at each\n", "step, both updates use the *current* $P(t)$ and $M(t)$ (not the\n", "newly-updated value of one when computing the other)." ], "id": "b5b38df4-ab99-4a0b-b858-9fb02820aa1b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 2.2 (13)\n", "\n", "Implement your update rules from Problem 2.1 and use them to simulate\n", "the lake under a loading of\n", "$L = 10\\ \\mu\\text{g}/(\\text{L}\\cdot\\text{yr})$, starting from\n", "$P(0) = M(0) = 0$.\n", "\n", "Two different quantities have to be “large enough” here, for two\n", "different reasons:\n", "\n", "- the **simulation horizon** $T$ has to be long enough for the slow\n", " sediment stock $M$ to fill. Its characteristic timescale is\n", " $1/b = 1/0.02 = 50\\ \\text{yr}$, so use $T = 300\\ \\text{yr}$ — six of\n", " those timescales. You will reuse this same run in Problem 2.3.\n", "- the **step size** $\\Delta t$ is what you are actually checking for\n", " numerical convergence, at that fixed $T$.\n", "\n", "Run a convergence study on $P(30)$, the phosphorus concentration thirty\n", "years in. Build a reference solution at a very fine step size (smaller\n", "than you think you might need), then test a sequence of\n", "successively-halved, coarser step sizes against it. Here is a starting\n", "point:\n", "\n", "``` julia\n", "function simulate_lake(loading, horizon, step)\n", " # TODO: forward Euler on P and M, using your update rules from Problem 2.1.\n", " # Both updates use the *current* P and M, not a partly-updated value.\n", " # Return the full P trajectory as a vector, one entry per step taken.\n", " error(\"simulate_lake has not been written yet -- replace this line with your loop\")\n", "end\n", "\n", "step_reference = 0.001 # 300,000 steps at T = 300; a few seconds\n", "P_reference = simulate_lake(10.0, 300.0, step_reference)\n", "P30_reference = P_reference[Int(30 / step_reference)] # the entry at t = 30 yr\n", "```\n", "\n", "Report a table of your step sizes, the resulting $P(30)$, the error\n", "relative to your reference, and the empirical order of convergence\n", "between successive refinements. Make a log-log plot of error\n", "vs. $\\Delta t$, and state which step size you will use for the rest of\n", "this problem and why.\n", "\n", "> **A step size that is too large will not just be inaccurate**\n", ">\n", "> Forward Euler can blow up entirely, rather than merely losing\n", "> accuracy, when the step is too large for the *fastest* process in the\n", "> system. Here that is $P$ itself, with a timescale of\n", "> $1/(s_\\text{out} + h) \\approx 1.7\\ \\text{yr}$ — much faster than the\n", "> sediment pool that set $T$. **The horizon is governed by the slowest\n", "> process; the step size is governed by the fastest one.**\n", ">\n", "> Start your sequence at $\\Delta t = 2\\ \\text{yr}$ and halve from there.\n", "> If you try something much coarser you may see the simulation diverge\n", "> or oscillate wildly; that is worth a sentence, but it is not what the\n", "> convergence study is measuring." ], "id": "5a72d587-f2ab-41fb-8f5a-b4c1e5a2f5f5" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 2.3 (10)\n", "\n", "Using the step size you justified in Problem 2.2, plot $P(t)$ over the\n", "full 300-year simulation.\n", "\n", "You already have $P(30)$ from Problem 2.2. What does $P$ look like once\n", "the simulation has fully run? What is happening physically to the lake\n", "between those two points in time, and around when does it happen?" ], "id": "fffe8019-6f1a-4253-9a1a-1ae03e6bfcf2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 2.4 (7)\n", "\n", "Suppose a lake manager instead used a simpler one-stock model that\n", "ignores the sediment pool entirely (as if $\\kappa = 0$, so there is no\n", "recycling at all):\n", "\n", "$$\\frac{dP}{dt} = L - s_\\text{out}P - hP.$$\n", "\n", "What is this simplified model’s equilibrium $P$ at\n", "$L = 10\\ \\mu\\text{g}/(\\text{L}\\cdot\\text{yr})$, and would the manager\n", "have concluded the lake is safely oligotrophic at this loading (recall\n", "that sediment release becomes significant once $P$ approaches\n", "$m = 25\\ \\mu\\text{g}/\\text{L}$)?\n", "\n", "Compare this to what you found in Problem 2.3. What does this tell you\n", "about the risk of using a simplified model, or a short simulation, to\n", "set a “safe” loading level? What would you recommend a lake manager do\n", "differently, given what you now know about this system?" ], "id": "42c85e19-0b7e-4e6b-8cf8-404dbfce57d9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "List any external references consulted, including classmates." ], "id": "c38afafb-3812-4cc3-819d-439e17226d8b" } ], "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" } } }