{ "cells": [ { "cell_type": "markdown", "id": "3d0a47e3-7c72-43fc-8563-d44cf886a24b", "metadata": {}, "source": [ "---\n", "title: \"A6: Transport Equation\"\n", "subtitle: \"Assignment 6: Transport Equation and the Courant–Friedrichs–Lewy (CFL) condition\"\n", "date: 2026-04-20\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: A6.ipynb\n", "---\n", "\n", ":::{.callout-note}\n", "\n", "* Complete the following and submit to Canvas before April 29, 23:59\n", "* Late work will recieve 0%,\n", "* Each assignment is worth the same, \n", "* Please get in contact with the TA/instructor in plenty of time if you need help,\n", "* Before submitting your work, make sure to check everything runs as expected. Click **Kernel > Restart Kernel and Run All Cells**.\n", "* Feel free to add more cells to experiment or test your answers,\n", "* I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own,\n", "* The use of GenAI is prohibited as outlined in the course syllabus. If I suspect you of cheating, you may be asked to complete a written or oral exam on the content of this assignment.\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 2, "id": "1f4163eb", "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra, Plots, LaTeXStrings" ] }, { "cell_type": "markdown", "id": "7f848e4a", "metadata": {}, "source": [ "::: {.box}\n", "\n", "***Name:***\n", "\n", "***Approx. time spent on this assignment:***\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "790bd5b2", "metadata": {}, "source": [ "In lectures, we have looked at \n", "\n", "\\begin{align}\n", " \\tag{IVP} &u'(t) = f(t, u(t)) \\nonumber\\\\\n", " \\tag{BVP} &-u'' + p u' + qu = f\\nonumber\\\\\n", " \\tag{heat} &u_t - \\alpha^2 \\nabla u = f\n", "\\end{align}\n", "\n", "and various methods to approximate the solution (IVP: Euler, one-step, Runge-Kutta, multistep; BVP: finite differences (including with FFT), finite elements; Heat equation: finite differences). " ] }, { "cell_type": "markdown", "id": "d980a92e", "metadata": {}, "source": [ "In this assignment, you'll analyse simple finite difference approximations to the linear transport equation:\n", "\n", "\\begin{align}\n", " &u_t + c u_x = 0 \\nonumber\\\\\n", " &u(x,0) = u_0(x) \\tag{transport}\n", "\\end{align}\n", "\n", "with constant velocity $c \\in \\mathbb R$. In lectures, we saw that the solution is given by $u(x,t) = u_0(x-ct)$. While in this case the exact solution is known, it is still useful to analyse methods for approximating this solution numerically. These methods apply more generally to more complicated hyperbolic equations. \n", "\n", "Consider the spatial domain $[a,b]$ and temporal domain $[0,T]$ with meshes given as follows: $x_i := a + \\frac{b-a}{m} i$ for $i = 0,\\dots,m$ and $t_j := j\\frac{T}{n}$ for $j = 0,\\dots,n$. Again, we define $h := \\frac{b-a}{m}$ and $\\tau := \\frac{T}{n}$. We aim to approximate $u(x_i, t_j)$ with $u_{i,j}$. \n", "\n", "**Exercise 1.** Suppose one approximates the linear transport equation with simple forward differences for both the space and time derivatives. Show that this approximation gives \n", "\n", "\\begin{align}\n", " u_{i,j+1} = (1+\\sigma) u_{i,j} - \\sigma u_{i+1,j}\n", " \\tag{transport$_{\\tau, h}$}\n", "\\end{align}\n", "\n", "where $\\sigma:= \\frac{c \\tau}{h}$." ] }, { "cell_type": "markdown", "id": "cca1d6d1", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "90ced6df", "metadata": {}, "source": [ "**Exercise 2.** We consider the linear transport equation with periodic boundary conditions: $u(a,t) = u(b,t)$ and $u_x(a,t) = u_x(b,t)$. Complete the following code sample to implement the finite difference scheme from exercise 1:" ] }, { "cell_type": "code", "execution_count": null, "id": "7caa8681", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "σ = -1.0\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\math5\\Math 5485\\Math5486\\pics\\transport-1.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\math5\\\\Math 5485\\\\Math5486\\\\pics\\\\transport-1.gif\")" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function fdm_transport(m, n, c, u0; a=0., b=1., T=5)\n", " h, τ = (b-a)/m, T/n \n", " x, t = collect(a:h:b), collect(0:τ:T)\n", " \n", " u = zeros(m+1, n+1)\n", " # Initial condition\n", " u[:, 1] = u0.(x)\n", "\n", " # CLF\n", " σ = c * τ / h \n", " println(\"σ = $(round(σ, digits=2))\")\n", "\n", " for j ∈ 1:n \n", " for i ∈ 1:m \n", " # periodic neighbour of m is 1:\n", " i′ = (i%m) + 1\n", "\n", " # finite difference rule here:\n", " u[i,j+1] = # your code here\n", " end\n", "\n", " # periodic boundary\n", " u[m+1,j+1] = # your code here\n", " end\n", "\n", " return x, t, u\n", "end\n", "\n", "a, b, T = -10, 10, 20\n", "m, n, c = 100, 100, -1.\n", "\n", "# initial condition (w/ PBCs)\n", "g(x) = sin(x) * exp(-x^2)\n", "u0(x) = g(a + mod(x - a, b-a))\n", "\n", "# exact solution \n", "u_e(x, t) = u0(x - c*t)\n", "\n", "x,t,u = fdm_transport(m, n, c, u0; a, b, T)\n", "\n", "anim = @animate for j ∈ 1:length(t)\n", " scatter(x, u[:, j]; \n", " title = L\"Transport equation $u_t + c u_{x} = 0$\",\n", " label = \"t = $(t[j])\", legend=:topright,\n", " ylims = [-1/2,1/2], m=2)\n", " plot!( x->u_e(x,t[j]); \n", " label=\"exact solution\", linestyle=:dot)\n", "end\n", "gif(anim, \"pics/transport-1.gif\", fps=10)\n" ] }, { "cell_type": "markdown", "id": "b5baab72", "metadata": {}, "source": [ "**Exercise 3.** Experiment with the above code. What happens when \n", "\n", "* *(i)* $\\sigma < -1$, \n", "* *(ii)* $\\sigma = -1$, \n", "* *(iii)* $-1 < \\sigma < 0$, \n", "* *(iv)* $\\sigma = 0$, and \n", "* *(v)* $\\sigma>0$.\n", "\n", "Either plot the numerical solution in each case, or write down the parameters you used." ] }, { "cell_type": "markdown", "id": "243f1580", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "4d72558b", "metadata": {}, "source": [ "The *domain of dependence* of the solution at $(x,t)$ is the set of all $y$ for which changing $u_0(y)$ can affect the solution $u$ at $(x,t)$. \n", "\n", "**Exercise 4.** Write down the domain of dependence of the solution to the linear transport equation at $(x,t)$." ] }, { "cell_type": "markdown", "id": "17ca13b9", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "fcfb241f", "metadata": {}, "source": [ "The *numerical domain of dependence* at $(x_i,t_j)$ is the set of all $\\{x_k\\}$ such that the initial data $u_{k,0}$ can affect $u_{i,j}$.\n", "\n", "**Exercise 5.** Write down the numerical domain of dependence of $\\text{(transport}_{h,\\tau}\\text{)}$." ] }, { "cell_type": "markdown", "id": "2170c8b0", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "5cb241e4", "metadata": {}, "source": [ "The *Courant–Friedrichs–Lewy (CFL) condition* states that, for a numerical method to converge, the exact domain of dependence must be contained in the numerical domain of dependence as $h,\\tau \\to 0$. \n", "\n", "**Exercise 6.** Under what condition on $\\sigma$ does our numerical method for the linear transport equation satisfy the CLF condition. " ] }, { "cell_type": "markdown", "id": "825f3cac", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "ea38bd22", "metadata": {}, "source": [ "**Remark:** If the CLF condition is violated, then one can change some of the initial data that is both outside the domain of numerical dependence (so **does not** affect the numerical solution) and inside the domain of exact dependence (so **does** affect the exact solution). As a result, the numerical scheme cannot approximate the true solution. The CLF condition is necessary but not sufficient for a numerical scheme to converge (as we will now see):" ] }, { "cell_type": "markdown", "id": "f6c22aa0", "metadata": {}, "source": [ "Consider the following centered difference approximation to the linear transport equation:\n", "\n", "\\begin{align}\n", " u_{i,j+1} = -\\frac12 \\sigma u_{i+1,j} + u_{ij} + \\frac{1}2\\sigma u_{i-1,j}\n", "\\end{align}\n", "\n", "where $\\sigma = \\frac{c\\tau}{h}$.\n", "\n", "**Exercise 7.** Show that, for $|\\sigma| \\leq 1$, this scheme satisfies the CLF condition." ] }, { "cell_type": "markdown", "id": "0f782853", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "6a778e99", "metadata": {}, "source": [ "**Exercise 8.** Show that this scheme is unstable for all $\\sigma\\not=0$ via a von Neumann stability analysis." ] }, { "cell_type": "markdown", "id": "385e55da", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**You answer here...**\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "1291576c", "metadata": {}, "source": [ "Here, we implement this method (and we can see that it is unstable for our choice of parameters): " ] }, { "cell_type": "code", "execution_count": 38, "id": "c6f68405", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "σ = 0.5\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\math5\\Math 5485\\Math5486\\pics\\transport-2.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\math5\\\\Math 5485\\\\Math5486\\\\pics\\\\transport-2.gif\")" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#| code-fold: true\n", "\n", "function fdm_transport_unstable(m, n, c, u0; a=0., b=1., T=5)\n", " h, τ = (b-a)/m, T/n \n", " x, t = collect(a:h:b), collect(0:τ:T)\n", " \n", " u = zeros(m+1, n+1)\n", " u[:, 1] = u0.(x)\n", "\n", " # CLF\n", " σ = c * τ / h \n", " println(\"σ = $(round(σ, digits=2))\")\n", "\n", " for j ∈ 1:n \n", " for i ∈ 1:m \n", " i₋ = (i==1) ? m : i-1\n", " i₊ = (i==m) ? 1 : i+1\n", " u[i,j+1] = (1/2)*σ*u[i₋,j] + u[i,j] + (-1/2)*σ*u[i₊,j] \n", " end\n", " u[m+1,j+1] = u[1,j+1]\n", " end\n", "\n", " return x, t, u\n", "end\n", "\n", "a, b, T = -10, 10, 20\n", "m, n, c = 100, 100, .5\n", "\n", "# initial condition (w/ PBCs)\n", "g(x) = sin(x) * exp(-x^2)\n", "u0(x) = g(a + mod(x - a, b-a))\n", "\n", "# exact solution \n", "u_e(x, t) = u0(x - c*t)\n", "\n", "x,t,u = fdm_transport_unstable(m, n, c, u0; a, b, T)\n", "\n", "anim = @animate for j ∈ 1:length(t)\n", " plot(x, u[:, j]; \n", " title = L\"Transport equation $u_t + c u_{x} = 0$ (unstable)\",\n", " label = \"t = $(t[j])\", legend=:topright,\n", " ylims = [-1/2,1/2], m=2)\n", " plot!( x->u_e(x,t[j]); \n", " label=\"exact solution\", linestyle=:dot)\n", "end\n", "gif(anim, \"pics/transport-2.gif\", fps=10)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.12", "language": "julia", "name": "julia-1.12" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }