{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab 1 Solutions\n", "\n", "> **Note**\n", ">\n", "> Labs are graded 0–2 on effort toward demonstrated work, not on\n", "> correctness. The notes in this key describe what a complete answer\n", "> might look like.\n", "\n", "## Setup\n", "\n", "The model setup is as given in the lab." ], "id": "35947c39-d790-49aa-bbb4-50cc8b328b56" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "landfill_params = let\n", " k = 0.15 # generation-response rate [1/yr]\n", " L0 = 100.0 # methane generation potential [m³ CH₄/Mg]\n", " Wmax = 1.0e5 # waste acceptance rate while open [Mg/yr]\n", " active_life = 20.0 # years the landfill accepts waste [yr]\n", " (; k, L0, Wmax, active_life)\n", "end\n", "\n", "T_total = 40.0 # total simulation horizon [yr]\n", "G0 = 0.0 # initial gas generation rate [m³ CH₄/yr]\n", "\n", "# Waste acceptance rate: constant while the landfill is open, zero after closure.\n", "function waste_rate(t, p)\n", " if t < p.active_life\n", " return p.Wmax\n", " else\n", " return 0.0\n", " end\n", "end" ], "id": "4" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 1: Discretize and Implement\n", "\n", "### Derivation\n", "\n", "Approximate the derivative at time $t$ with a forward difference over\n", "one step:\n", "\n", "$$\\frac{dG}{dt} \\approx \\frac{G(t + \\Delta t) - G(t)}{\\Delta t}.$$\n", "\n", "Forward Euler evaluates the right-hand side at the *start* of the step,\n", "using what we already know at time $t$:\n", "\n", "$$\\frac{G(t + \\Delta t) - G(t)}{\\Delta t} = k\\left(L_0 W(t) - G(t)\\right).$$\n", "\n", "Solving for the unknown:\n", "\n", "$$G(t + \\Delta t) = G(t) + \\Delta t \\, k\\left(L_0 W(t) - G(t)\\right).$$\n", "\n", "We can consolidate:\n", "\n", "$$G(t + \\Delta t) = \\left(1 - k\\Delta t\\right) G(t) + k\\Delta t \\, L_0 W(t).$$\n", "\n", "Each step blends the current gas generation rate with the rate the waste\n", "in place would eventually produce. $L_0 W(t)$, and $k\\Delta t$ determine\n", "the influence of new waste.\n", "\n", "> **Note**\n", ">\n", "> Evaluating $W$ or $G$ at $t + \\Delta t$ on the right-hand side is\n", "> wrong (that is no longer forward Euler). For forward Euler, the\n", "> **present state** (at time $t$) is what determines the update rule.\n", "\n", "### Implementation\n", "\n", "The only line to fill in is the update, which is the derivation above\n", "written in code." ], "id": "131cea63-7d19-430f-9580-ac9a59cf455d" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "landfill_gas_rate(G, t, p) = p.k * (p.L0 * waste_rate(t, p) - G)\n", "\n", "function landfill_gas_simulate(G_ic, T, Δt, p)\n", " steps = Int(round(T / Δt))\n", " G = zeros(steps + 1) # index 1 holds the initial condition\n", " G[1] = G_ic\n", " for i in 1:steps\n", " t = (i - 1) * Δt\n", " G[i+1] = G[i] + Δt * landfill_gas_rate(G[i], t, p)\n", " end\n", " return G\n", "end" ], "id": "6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sanity check:" ], "id": "74f91747-9632-4d49-a244-947ff3d969a7" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "Δt_check = 0.5\n", "G_check = landfill_gas_simulate(G0, T_total, Δt_check, landfill_params)\n", "times_check = collect(0:length(G_check)-1) .* Δt_check\n", "\n", "p_check = plot(times_check, G_check ./ 1e6, color=cb_blue, legend=false,\n", " xlabel=\"Time [yr]\", ylabel=\"Gas generation [million m³/yr]\")\n", "vline!(p_check, [landfill_params.active_life], color=:black, linewidth=2,\n", " linestyle=:dash)\n", "plot!(p_check, size=(900, 380), left_margin=10mm, bottom_margin=10mm)" ], "id": "cell-fig-sanity" }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shape makes sense given what the model tells us. While the landfill\n", "is open (from $t=0$ to $t=20$), $G$ rises quickly at first and then\n", "flattens as it approaches $L_0 W_\\text{max} =$ 10,000,000 m3\n", "CH4/yr, the rate that constant waste acceptance would\n", "eventually support. It never gets there: by facility closure it has\n", "reached 96% of that level. The peak falls exactly at closure, because\n", "that is the moment the source switches off. After closure, $G$ decays\n", "exponentially on the same timescale it rose, down to 4% of the peak by\n", "year 40.\n", "\n", "## Problem 2: How Small Does $\\Delta t$ Need to Be?\n", "\n", "### Problem 2a: Characteristic Timescale\n", "\n", "The model has one rate, $k = 0.15\\ \\text{yr}^{-1}$, so the\n", "characteristic timescale is\n", "\n", "$$\\frac{1}{k} = \\frac{1}{0.15\\ \\text{yr}^{-1}} \\approx 6.7\\ \\text{yr}$$\n", "\n", "The waste acceptance rate $W(t)$ is a forcing, not a rate. It switches\n", "off abruptly at closure, but $G$ still responds to that switch on the\n", "$1/k$ timescale.\n", "\n", "For a starting step, we might pick $\\Delta t = 0.5$ yr. That is about\n", "thirteen steps per timescale, which is well inside it, and it divides 20\n", "evenly, so a grid point lands exactly on closure. But any value\n", "reasonably below the timescale is defensible as a starting point.\n", "\n", "> **Impact of choosing too large an initial step**\n", ">\n", "> Look at the rearranged update rule from Problem 1. The weight on\n", "> $G(t)$ is $1 - k\\Delta t$. Once $\\Delta t$ exceeds $1/k$, that weight\n", "> is negative: after closure, each step flips the sign of $G$, and the\n", "> model reports *negative* gas generation, which is physically\n", "> impossible. Beyond $2/k \\approx 13$ yr the flips grow and the solution\n", "> blows up. So $1/k$ is a hard ceiling on $\\Delta t$, and you usually\n", "> want to pick something much smaller than that." ], "id": "8a161a5a-0c16-4741-b8a9-3897ca1ede58" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "too_coarse = landfill_gas_simulate(G0, T_total, 7.0, landfill_params)\n", "most_negative = minimum(too_coarse)" ], "id": "12" }, { "cell_type": "markdown", "metadata": {}, "source": [ "At $\\Delta t = 7$ yr, just past the ceiling, the lowest value of $G$ is\n", "−500,062 m³ CH4/yr.\n", "\n", "### Problem 2b: Convergence Results\n", "\n", "Let’s build the reference. A step of $10^{-4}$ yr is 400,000 steps,\n", "which still runs in a few milliseconds for this simple model." ], "id": "f0c87b54-a796-40d2-bc7b-96ae4bab384d" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "Δt_ref = 1e-4\n", "G_ref = landfill_gas_simulate(G0, T_total, Δt_ref, landfill_params)\n", "peak_ref = maximum(G_ref)" ], "id": "16" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The reference peak is 9,502,141 m3 CH4/yr.\n", "\n", "> **Checking the Reference**\n", ">\n", "> You did not have to do this, but here we can actually find an analytic\n", "> solution to the maximum gas generation to compare the reference to\n", "> (this is not possible in general). While the landfill is open,\n", "> $G(t) = L_0 W_\\text{max}\\left(1 - e^{-kt}\\right)$, so the true peak at\n", "> closure is\n", ">\n", "> $$L_0 W_\\text{max}\\left(1 - e^{-k \\cdot 20}\\right) = 10^7 \\left(1 - e^{-3}\\right),$$\n", ">\n", "> which is 9,502,129 m³ CH4/yr.\n", ">\n", "> The reference is off by 11 m³ CH4/yr, far smaller than any\n", "> error in the table below, so it is safe to measure against. Again, if\n", "> you don’t have an analytic solution, just pick something that seems\n", "> very small.\n", "\n", "Starting from the Problem 2a estimate and halving four times:" ], "id": "683635d9-af89-4484-95be-4a410b7cbea1" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "Δts = [0.5, 0.25, 0.125, 0.0625, 0.03125]\n", "peaks = zeros(length(Δts))\n", "errs = zeros(length(Δts))\n", "for i in 1:length(Δts)\n", " G_test = landfill_gas_simulate(G0, T_total, Δts[i], landfill_params)\n", " peaks[i] = maximum(G_test)\n", " errs[i] = abs(peaks[i] - peak_ref)\n", "end\n", "\n", "orders = zeros(length(Δts) - 1)\n", "for i in 1:length(Δts)-1\n", " orders[i] = log2(errs[i] / errs[i+1])\n", "end" ], "id": "18" }, { "cell_type": "markdown", "metadata": {}, "source": [ "| $\\Delta t$ \\[yr\\] | Peak $G$ \\[m³ CH₄/yr\\] | Error vs. reference | Error improvement |\n", "|-----------------:|-----------------:|-----------------:|-----------------:|\n", "| 0.5 | 9,557,749 | 55,608 | — |\n", "| 0.25 | 9,530,042 | 27,901 | 0.99 |\n", "| 0.125 | 9,516,109 | 13,969 | 1.00 |\n", "| 0.0625 | 9,509,125 | 6,985 | 1.00 |\n", "| 0.03125 | 9,505,629 | 3,488 | 1.00 |" ], "id": "443345b7-3ba8-4c33-b9fb-a3e732df6e1f" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "error_ticks = [5e3, 1e4, 2e4, 5e4]\n", "p_conv = plot(Δts, errs, xscale=:log10, yscale=:log10, markershape=:circle,\n", " markersize=7, color=cb_vermillion, label=\"forward Euler\",\n", " xticks=(Δts, string.(Δts)), yticks=(error_ticks, with_commas.(error_ticks)),\n", " xlabel=\"Δt [yr]\", ylabel=\"Error in peak [m³/yr]\", legend=:bottomright)\n", "# Slope-1 reference: errors parallel to it mean a first-order method.\n", "plot!(p_conv, Δts, errs[1] .* (Δts ./ Δts[1]), color=:black, linewidth=2,\n", " linestyle=:dash, label=\"slope 1\")\n", "plot!(p_conv, size=(900, 400), left_margin=12mm, bottom_margin=10mm)" ], "id": "cell-fig-convergence" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The error halves every time $\\Delta t$ halves: each improvement is close\n", "to 1, and the errors run parallel to the slope-1 line. That is\n", "first-order convergence, which is what the Taylor argument predicts.\n", "Forward Euler keeps the first two terms of\n", "$G(t + \\Delta t) = G(t) + \\Delta t \\, G'(t) + \\frac{1}{2}\\Delta t^2 G''(t) + \\cdots$,\n", "so each step commits an error proportional to $\\Delta t^2$. Reaching\n", "closure takes $20/\\Delta t$ steps, and $20/\\Delta t$ errors of size\n", "$\\Delta t^2$ add up to an error proportional to $\\Delta t$.\n", "\n", "The errors also all have the same sign: forward Euler *over*-predicts\n", "the peak at every step size. While the landfill fills, $dG/dt$ shrinks\n", "over each step as $G$ climbs toward $L_0 W_\\text{max}$. Forward Euler\n", "uses the rate from the start of the step, when it is largest, so it\n", "overshoots a little on every step.\n", "\n", "> **Some Pathologies Unrelated to Coding**\n", ">\n", "> You might see the following, but these are not code errors:\n", ">\n", "> - **Improvements that wobble around 1.** If you start from $1/(10k)$\n", "> and rounded it to 0.67 yr get a sequence whose grid points never\n", "> land on closure. The step that straddles year 20 keeps adding waste\n", "> for part of a step, and that fraction changes irregularly as\n", "> $\\Delta t$ halves. Starting from 0.67, the improvements run 0.87,\n", "> 0.77, 1.30, 0.84, while the slope across the whole sequence is 0.95.\n", "> This is still first order.\n", "> - **Improvements that creep above 1 at the fine end.** This one *is* a\n", "> problem, with the reference rather than the method:\n", "> $\\Delta t_\\text{ref}$ was not much finer than the finest test step.\n", "> The reference carries an error of the same sign, so the measured\n", "> errors shrink faster than the true ones. With\n", "> $\\Delta t_\\text{ref} = 0.01$ yr and the step sizes above, the\n", "> improvements run 1.02, 1.06, 1.13, 1.30. The fix is a finer\n", "> reference.\n", "\n", "## Problem 3: Make and Justify a Recommendation\n", "\n", "From this analysis, we might recommend using $\\Delta t = 0.25$ yr, a\n", "quarterly step.\n", "\n", "Other step sizes are equally defensible. Even an annual step puts the\n", "peak only 1.2% high, which is a real feature of this model: a timescale\n", "of nearly seven years is forgiving." ], "id": "283203d2-2e04-4d4f-8737-8d2715908c42" } ], "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" } } }