{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Partial differential equations\n", "\n", "Partial differential equations occur when we are dealing with functions of more than one variable, for instance *fields*.\n", "\n", "Examples:\n", "\n", "- Electrostatic potential $\\phi (x,y,z)$ (Poisson's equation)\n", "$$\n", "\\Delta\\phi(x,y,z) = -\\frac{\\rho(x,y,z)}{\\epsilon_0}\n", "$$\n", "\n", "- Density or temperature profiles (diffusion/heat equation)\n", "$$\n", "\\frac{\\partial u(\\mathbf{x},t)}{\\partial t} = D \\Delta u(\\mathbf{x},t)\n", "$$\n", "\n", "- Displacement (amplitude) profile (wave equation)\n", "$$\n", "\\frac{\\partial^2 u(\\mathbf{x},t)}{\\partial t^2} = c^2\\Delta u(\\mathbf{x},t)\n", "$$\n", "\n", "- Fluid dynamical fields (flow velocity) -- e.g. Navier-Stokes equations " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## General methods for solving PDE's\n", "\n", "- **Finite difference method**\n", " - Approximate the derivatives by finite differences\n", " - Easier to implement than other methods\n", " - Works best for regular (rectangular) shapes\n", "- Finite element method\n", " - Subdivide the system into smaller parts -- finite elements\n", " - Boundary value problems in 2/3 dimensions\n", " - Works well for irregular shapes\n", "- Finite volume method\n", " - Convert surface integrals around each mesh point into volume integrals\n", " - Conserves mass by design, good for fluid dynamical equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary value problems\n", "\n", "\"Electrostatics \n", "
\n", " Source: M. Newman, Computational Physics\n", "
\n", "\n", "Boundary value problems typically deal with static solutions (no time variable)\n", "\n", "Consider Laplace's equation (no external charges) in two dimensions\n", "$$\n", "\\frac{\\partial^2 \\phi(x,y)}{\\partial x^2}\n", "+\n", "\\frac{\\partial^2 \\phi(x,y)}{\\partial y^2}\n", "=\n", "0.\n", "$$\n", "\n", "\n", "\n", "Boundary conditions: one of the walls has voltage $V$ applied to it, the rest are grounded:\n", "\\begin{align*}\n", "\\phi(x,L) & = V, \\\\\n", "\\phi(x,0) & = 0, \\\\\n", "\\phi(0,y) & = 0, \\\\\n", "\\phi(L,y) & = 0.\n", "\\end{align*}\n", "\n", "Let us apply finite difference method. \n", "Discretize the box into a grid in steps $a = L / M$ in each direction.\n", "Approximate the derivatives by central differences:\n", "\\begin{align*}\n", "\\frac{\\partial^2 \\phi(x,y)}{\\partial x^2} & = \\frac{\\phi(x+a,y) - 2\\phi(x,y) + \\phi(x-a,y)}{a^2}, \\\\\n", "\\frac{\\partial^2 \\phi(x,y)}{\\partial y^2} & = \\frac{\\phi(x,y+a) - 2\\phi(x,y) + \\phi(x,y-a)}{a^2}.\n", "\\end{align*}\n", "\n", "The Laplace's equation becomes\n", "$$\n", "\\phi(x+a,y) + \\phi(x-a,y) + \\phi(x,y+a) + \\phi(x,y-a) - 4 \\phi(x,y) = 0.\n", "$$\n", "\n", "In principle, this is a system of $N \\sim M^2$ linear equations that can be solved exactly.\n", "The general method will however require $O(N^3) = O(M^6)$ operations which becomes unfeasible already for $M \\sim 100$ or so." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Jacobi (relaxation) method\n", "\n", "Start with some initial guess $\\phi_0(x,y)$. The next iteration is given by\n", "$$\n", "\\phi_{n+1}(x,y) = \\frac{\\phi_n(x+a,y) + \\phi_n(x-a,y) + \\phi_n(x,y+a) + \\phi_n(x,y-a)}{4}\n", "$$\n", "for all points inside the box. Preserve the boundary conditions at each iteration:\n", "\\begin{align*}\n", "\\phi_{n+1}(x,L) & = V, \\\\\n", "\\phi_{n+1}(x,0) & = 0, \\\\\n", "\\phi_{n+1}(0,y) & = 0, \\\\\n", "\\phi_{n+1}(L,y) & = 0.\n", "\\end{align*}\n", "\n", "This is the Jacobi (or relaxation) method, and for Laplace's equation it is always converges." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the Jacobi method\n", "# The new field is written into phinew\n", "def iteration_jacobi(phinew, phi):\n", " M = len(phi) - 1\n", " \n", " # Boundary conditions\n", " phinew[0,:] = phi[0,:]\n", " phinew[M,:] = phi[M,:]\n", " phinew[:,0] = phi[:,0]\n", " phinew[:,M] = phi[:,M]\n", " \n", " for i in range(1,M):\n", " for j in range(1,M):\n", " phinew[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4\n", " \n", " delta = np.max(abs(phi-phinew))\n", " \n", " return delta\n", "\n", "def jacobi_solve(phi0, target_accuracy = 1e-6, max_iterations = 100):\n", " delta = target_accuracy + 1.\n", " phi = phi0.copy()\n", " for i in range(max_iterations):\n", " delta = iteration_jacobi(phi, phi0)\n", " phi0, phi = phi, phi0\n", " \n", " if (delta <= target_accuracy):\n", " print(\"Jacobi method converged in \" + str(i+1) + \" iterations\")\n", " return phi0\n", " \n", " print(\"Jacobi method failed to converge to a required precision in \" + str(max_iterations) + \" iterations\")\n", " print(\"The error estimate is \", delta)\n", " \n", " return phi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# Measure execution time\n", "\n", "# Constants\n", "M = 100 # Grid squares on a side\n", "V = 1.0 # Voltage at top wall\n", "target = 1e-4 # Target accuracy\n", "\n", "# Initialize with zeros\n", "phi = np.zeros([M+1,M+1],float)\n", "# Boundary condition\n", "phi[0,:] = V\n", "phi[:,0] = 0\n", "\n", "phi = jacobi_solve(phi, target, 10000)\n", "\n", "# Plot\n", "import matplotlib.pyplot as plt\n", "\n", "plt.title(\"Electrostatic potential\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "CS = plt.imshow(phi, vmax=1., vmin=0.,origin=\"upper\",extent=[0,M,0,M])\n", "plt.colorbar(CS)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "## This is adapted Example 9.1 from M. Newman, \"Computational Physics\"\n", "## Source: http://www-personal.umich.edu/~mejn/cp/programs/laplace.py\n", "\n", "from numpy import empty,zeros,max\n", "from pylab import imshow,gray,show\n", "\n", "# Constants\n", "M = 100 # Grid squares on a side\n", "V = 1.0 # Voltage at top wall\n", "target = 1e-4 # Target accuracy\n", "\n", "# Create arrays to hold potential values\n", "phi = zeros([M+1,M+1],float)\n", "phi[0,:] = V\n", "phiprime = empty([M+1,M+1],float)\n", "\n", "# Main loop\n", "delta = 1.0\n", "while delta>target:\n", "\n", " # Calculate new values of the potential\n", " for i in range(M+1):\n", " for j in range(M+1):\n", " if i==0 or i==M or j==0 or j==M:\n", " phiprime[i,j] = phi[i,j]\n", " else:\n", " phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \\\n", " + phi[i,j+1] + phi[i,j-1])/4\n", "\n", " # Calculate maximum difference from old values\n", " delta = max(abs(phi-phiprime))\n", "\n", " # Swap the two arrays around\n", " phi,phiprime = phiprime,phi\n", "\n", "# Make a plot\n", "# imshow(phi)\n", "# gray()\n", "# show()\n", "plt.title(\"Electrostatic potential\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "CS = plt.imshow(phi, vmax=1., vmin=0.,origin=\"upper\",extent=[0,M,0,M])\n", "plt.colorbar(CS)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss-Seidel method with overrelaxation\n", "\n", "The base Jacobi method corresponds to an iteration\n", "$$\n", "\\phi_{n+1}(x,y) = \\frac{\\phi_n(x+a,y) + \\phi_n(x-a,y) + \\phi_n(x,y+a) + \\phi_n(x,y-a)}{4}\n", "$$\n", "\n", "This requires to have two arrays independently.\n", "In Gauss-Seidel method one uses the already computed values of $\\phi_{n+1}$ where availabl instead.\n", "The method thus corresponds to\n", "$$\n", "\\phi_{n+1}(x,y) = \\frac{\\phi_n(x+a,y) + \\phi_{n+1}(x-a,y) + \\phi_n(x,y+a) + \\phi_{n+1}(x,y-a)}{4}\n", "$$\n", "\n", "Another modification is the use of overrelaxation.\n", "\n", "The base Jacobi method is a type of relaxation method\n", "$$\n", "\\phi_{n+1}(x,y) = \\phi_n(x,y) + \\Delta_n \\phi(x,y).\n", "$$\n", "\n", "A way to speed-up calculation is to overrelaxate the solution a bit\n", "$$\n", "\\phi_{n+1}(x,y) = \\phi_n(x,y) + (1+\\omega) \\Delta_n \\phi(x,y),\n", "$$\n", "where $\\omega > 0$.\n", "\n", "This implies the following iterative procedure\n", "$$\n", "\\phi_{n+1}(x,y) = (1+\\omega) \\frac{\\phi_n(x+a,y) + \\phi_n(x-a,y) + \\phi_n(x,y+a) + \\phi_n(x,y-a)}{4} - \\omega \\, \\phi_n(x,y).\n", "$$\n", "\n", "This is unstable for $\\omega > 0$. However, applied to Gauss-Seidel method, this corresponds to\n", "$$\n", "\\phi_{n+1}(x,y) = (1+\\omega) \\frac{\\phi_n(x+a,y) + \\phi_{n+1}(x-a,y) + \\phi_n(x,y+a) + \\phi_{n+1}(x,y-a)}{4} - \\omega \\, \\phi_n(x,y),\n", "$$\n", "and provides generally stable solution for $\\omega < 1$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the Jacobi method\n", "# The new field is written into phinew\n", "# omega >= 0 is the overelaxation parameter\n", "def gaussseidel_iteration(phi, omega = 0):\n", " M = len(phi) - 1\n", " \n", " delta = 0.\n", " \n", " # New iteration\n", " for i in range(1,M):\n", " for j in range(1,M):\n", " phiold = phi[i,j]\n", " phi[i,j] = (1. + omega) * (phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4 - omega * phi[i,j]\n", " delta = np.maximum(delta, abs(phiold - phi[i,j]))\n", " \n", " return delta\n", "\n", "def gaussseidel_solve(phi0, omega = 0, target_accuracy = 1e-6, max_iterations = 100):\n", " delta = target_accuracy + 1.\n", " phi = phi0.copy()\n", " for i in range(max_iterations):\n", " delta = gaussseidel_iteration(phi, omega)\n", " \n", " if (delta <= target_accuracy):\n", " print(\"Jacobi method converged in \" + str(i+1) + \" iterations\")\n", " return phi\n", " \n", " print(\"Jacobi method failed to converge to a required precision in \" + str(max_iterations) + \" iterations\")\n", " print(\"The error estimate is \", delta)\n", " \n", " return phi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "M = 100 # Grid squares on a side\n", "V = 1.0 # Voltage at top wall\n", "target = 1e-4 # Target accuracy\n", "\n", "# Initialize with zeros\n", "phi = np.zeros([M+1,M+1],float)\n", "# Boundary condition\n", "phi[0,:] = V\n", "phi[:,0] = 0\n", "\n", "\n", "omega = 0.93\n", "phi = gaussseidel_solve(phi, omega, target, 1000)\n", "\n", "# Plot\n", "plt.title(\"Electrostatic potential\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "CS = plt.imshow(phi, vmax=1., vmin=0.,origin=\"upper\",extent=[0,M,0,M])\n", "plt.colorbar(CS)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial value problem\n", "\n", "In most cases we are studying the time evolution of a certain profile $u(t,\\mathbf{x})$.\n", "In this case PDEs desribe the time evolution of the field(s) starting from some initial conditions and obeying boundary conditions.\n", "\n", "Let us take the heat equation as an example.\n", "The field $u$ is the temperature $T$.\n", "In one dimension the heat equation reads\n", "$$\n", "\\frac{\\partial u}{\\partial t} = D \\, \\frac{\\partial^2 u}{\\partial x^2},\n", "$$\n", "where $D$ is the thermal diffusivity constant.\n", "\n", "This equation describes the time evolution of $u(t,x)$ given initial profile\n", "$$\n", "u(t=0,x) = u_0(x),\n", "$$\n", "and boundary conditions\n", "\\begin{align*}\n", "u(t,x=0) & = u_{\\rm left}(t), \\\\\n", "u(t,x=L) & = u_{\\rm right}(t).\n", "\\end{align*}\n", "\n", "If $u_{\\rm left}(t)$ and $u_{\\rm right}(t)$ are time-independent, we know that the solution will approach a stationary profile as $t \\to \\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finite difference methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. FTCS scheme\n", "\n", "FTCS (Finite Time Centered Space) scheme is the simplest method for solving the heat equation\n", "\n", "First we discretize the spatial coordinate into a grid with $N + 1$ points, i.e.\n", "$$\n", "x_k = a k, \\qquad k = 0\\ldots N, \\qquad a = L/N,\n", "$$\n", "and approximate the derivative $\\partial^2 u / \\partial x^2$ by the lowest order central difference\n", "$$\n", "\\frac{\\partial^2 u(t,x)}{\\partial x^2} \\approx \\frac{u(t,x+a) - 2u(t,x) + u(t,x-a)}{a^2}.\n", "$$\n", "\n", "To evaluate the time evolution we will work is small time steps of size $h$.\n", "The time derivative is approximated by the forward difference\n", "$$\n", "\\frac{\\partial u(t,x)}{\\partial t} \\approx \\frac{u(t+h,x) - u(t,x)}{h}.\n", "$$\n", "\n", "This gives the following discretized equation\n", "$$\n", "\\frac{u(t+h,x) - u(t,x)}{h} = D \\frac{u(t,x+a) - 2u(t,x) + u(t,x-a)}{a^2}.\n", "$$\n", "The method is explicit: to evaluate $u(t+h,x)$ at the next time step we only need to know $u(t,x)$ profile at the present time step.\n", "Denoting the discretized time variable by superscript $n$ (such that $t_n = hn$) and the spatial variable by subscript $k$ (such that $x_k = ak$) we get the following iterative procedure\n", "$$\n", "u^{n+1}_k = u^n_k + r \\, (u^n_{k+1} - 2u^n_k + u^n_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "$$\n", "Here\n", "$$\n", "r \\equiv \\frac{Dh}{a^2}\n", "$$\n", "is a dimensionless parameter." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the FTCS scheme in the time direction\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def heat_FTCS_iteration(u, r):\n", " N = len(u) - 1\n", " \n", " unew = np.empty_like(u)\n", " \n", " # Boundary conditions\n", " unew[0] = u[0]\n", " unew[N] = u[N]\n", " \n", " \n", " # FTCS scheme\n", " for i in range(1,N):\n", " unew[i] = u[i] + r * (u[i+1] - 2 * u[i] + u[i-1])\n", " \n", " return unew\n", "\n", "\n", "# Perform nsteps FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def heat_FTCS_solve(u0, h, nsteps, a, D = 1.):\n", " u = u0.copy()\n", " r = h * D / a**2\n", " for i in range(nsteps):\n", " u = heat_FTCS_iteration(u, r)\n", " \n", " return u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us consider Example 9.3 from M. Newman, *Computational Physics*:\n", "\n", "We have a 1 cm long steel container, initially at a temperature 20° C.\n", "It is placed in bath of cold water at 0° C and filled on top with hot water at 50° C.\n", "Our goal is to calculate the temperature profile as function of time.\n", "The thermal diffusivity constant for stainless steel is $D = 4.25 \\cdot 10^{-6}$ m$^2$ s$^{-1}$.\n", "\n", "We will calculate the profile at times $t = 0.01$ s, $0.1$ s, $0.4$ s, $1$ s, and $10$ s." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-4 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with FTCS scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0, 10.0 ]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_FTCS_solve(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "plt.title(\"Heat equation\")\n", "plt.xlabel('${x}$ [cm]')\n", "plt.ylabel('${T}$ [Celsius]')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]) + \" s\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try a larger time step such that $r > 1/2$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1.2e-3 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with FTCS scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0, 10.0 ]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_FTCS_solve(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "plt.title(\"Heat equation\")\n", "plt.xlabel('${x}$ [cm]')\n", "plt.ylabel('${T}$ [Celsius]')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]) + \" s\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we increase the spatial size step to bring $r$ back below 1/2?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 90 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1.2e-3 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with FTCS scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0, 10.0 ]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_FTCS_solve(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "plt.title(\"Heat equation\")\n", "plt.xlabel('${x}$ [cm]')\n", "plt.ylabel('${T}$ [Celsius]')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]) + \" s\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Implicit scheme\n", "\n", "In the implicit scheme one uses backward difference for the time derivative.\n", "This implies\n", "$$\n", "\\frac{\\partial u(t+h,x)}{\\partial t} \\approx \\frac{u(t+h,x) - u(t,x)}{h},\n", "$$\n", "thus\n", "$$\n", "\\frac{u(t+h,x) - u(t,x)}{h} = D \\frac{u(t+h,x+a) - 2u(t+h,x) + u(t+h,x-a)}{a^2}.\n", "$$\n", "In discretized notation\n", "$$\n", "u^{n+1}_k = u^n_k + r \\, (u^{n+1}_{k+1} - 2u^{n+1}_k + u^{n+1}_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "$$\n", "In other words, we have a system of linear equations for $u^{n+1}_i$:\n", "$$\n", "-r u^{n+1}_{k-1} + (1+2r) u^{n+1}_{k} - r u^{n+1}_{k+1} = u^n_k, \\qquad k = 1 \\ldots N-1.\n", "$$\n", "The system should be solved at each time step. Luckily, the system is tridiagonal, so it is solved in linear time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Solve tridiagonal system of linear equations\n", "# d: vector of diagonal elements\n", "# l: vector of elements on the lower subdiagonal\n", "# u: vector of elements on the upper superdiagonal\n", "# v0: right-hand-side vector\n", "def linsolve_tridiagonal(d, l, u, v0):\n", " # Initialization\n", " N = len(v0)\n", " a = d.copy() # Current diagonal elements\n", " b = u.copy() # Current upper diagonal elements\n", " v = v0.copy()\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " if (a[r] == 0.):\n", " print(\"Diagonal element is zero! Cannot solve the tridiagonal system with simple Gaussian elimination\")\n", " return None\n", " b[r] /= a[r]\n", " v[r] /= a[r]\n", " a[r] = 1.\n", " if (r < N - 1):\n", " a[r + 1] -= l[r+1] * b[r]\n", " v[r + 1] -= l[r+1] * v[r]\n", " \n", " # Backsubstitution\n", " x = np.empty(N,float)\n", " \n", " x[N - 1] = v[N - 1]\n", " for r in range(N-2,-1,-1):\n", " x[r] = v[r] - b[r] * x[r + 1]\n", " \n", " return x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "\n", "# Single iteration of the FTCS scheme in the time direction\n", "# The new field is written into unew\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def heat_implicit_iteration(u, r):\n", " N = len(u) - 1\n", " \n", " unew = np.empty_like(u)\n", " \n", " # Boundary conditions\n", " unew[0] = u[0]\n", " unew[N] = u[N]\n", " \n", " d = np.full(N-1, 1+2.*r)\n", " ud = np.full(N-1, -r)\n", " ld = np.full(N-1, -r)\n", " v = np.array(u[1:N])\n", " v[0] += r * u[0]\n", " v[N-2] += r * u[N]\n", " \n", " unew[1:N] = linsolve_tridiagonal(d,ld,ud,v)\n", " \n", " return unew\n", "\n", "\n", "# Perform nsteps FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def heat_implicit_solve(u0, h, nsteps, a, D = 1.):\n", " u = u0.copy()\n", " r = h * D / a**2\n", " # print(\"Heat equation with r =\", r)\n", " for i in range(nsteps):\n", " u = heat_implicit_iteration(u, r)\n", " \n", " return u" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with implicit scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0, 10.0 ]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_implicit_solve(u, h, nsteps, a, D)\n", " # u = heat_FTCS_solve(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "plt.title(\"Heat equation\")\n", "plt.xlabel('${x}$ [cm]')\n", "plt.ylabel('${T}$ [Celsius]')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]) + \" s\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Animate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "from matplotlib.animation import FuncAnimation\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-3 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with implicit scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Heat equation\", fontsize = 18)\n", "current_time = 0\n", "def animate_heat_equation_1D(i):\n", " ax.clear()\n", " global fps\n", " niters = round((0.1/fps)/h)\n", " global u, current_time\n", " u = heat_implicit_solve(u, h, niters, a, D)\n", " current_time += h * niters\n", " # print(current_time)\n", " \n", " # plt.title(\"Heat equation\")\n", " plt.xlabel('${x}$ [cm]')\n", " plt.ylabel('${T}$ [Celsius]')\n", " plt.plot(xk,u,label=\"${t=}$\" + \"{:.3f}\".format(current_time) + \" s\")\n", " plt.legend()\n", " # plt.show()\n", "\n", "fps = 30\n", "tend = 2.\n", "ani = FuncAnimation(fig, animate_heat_equation_1D, frames=fps * round(tend) * 10, interval=1000/fps, repeat=False)\n", "#plt.show()\n", "\n", "ani.save(\"heat_equation_1D.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"heat_equation_1D.gif\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Crank-Nicolson scheme\n", "\n", "Crank-Nicolson method is a combination of FTCS and implicit schemes.\n", "Essentially it corresponds to approximating the time derivative as average of explicit and implicit methods\n", "$$\n", "\\frac{\\partial u(t,x)}{\\partial t} \\approx \\frac{1}{2} \\left[ D \\, \\frac{\\partial^2 u(t+h,x)}{\\partial x^2} + D \\, \\frac{\\partial^2 u(t,x)}{\\partial x^2} \\right].\n", "$$\n", "\n", "Applying central differences to the spatial derivatives we obtain\n", "$$\n", "\\frac{u(t+h,x) - u(t,x)}{h} = \\frac{D}{2} \\frac{u(t+h,x+a) - 2u(t+h,x) + u(t+h,x-a)}{a^2} + \\frac{D}{2} \\frac{u(t,x+a) - 2u(t,x) + u(t,x-a)}{a^2}.\n", "$$\n", "\n", "This corresponds to the following tri-diagonal system of linear equations\n", "$$\n", "-r u^{n+1}_{k-1} + 2(1+r) u^{n+1}_{k} - r u^{n+1}_{k+1} = ru^n_{k-1} + 2(1-r)u^n_k + ru^n_{k+1}, \\qquad k = 1 \\ldots N-1.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def heat_crank_nicolson_iteration(u, r):\n", " N = len(u) - 1\n", " \n", " unew = np.empty_like(u)\n", " \n", " # Boundary conditions\n", " unew[0] = u[0]\n", " unew[N] = u[N]\n", " \n", " d = np.full(N-1, 2*(1+r))\n", " ud = np.full(N-1, -r)\n", " ld = np.full(N-1, -r)\n", " \n", " # Crank-Nicolson explicit step\n", " v = u[1:N]*2*(1-r) + u[:-2]*r + u[2:]*r\n", " v[0] += r * u[0]\n", " v[N-2] += r * u[N]\n", " \n", " unew[1:N] = linsolve_tridiagonal(d, ld, ud, v)\n", " \n", " return unew\n", "\n", "def heat_crank_nicolson_solve(u0, h, nsteps, a, D = 1.):\n", " u = u0.copy()\n", " r = h * D / (a**2)\n", " for i in range(nsteps):\n", " u = heat_crank_nicolson_iteration(u, r)\n", " \n", " return u\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with Crank-Nicolson scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1],float)\n", "# Initial temperature\n", "u[1:N] = Tmid\n", "# Boundary conditions\n", "u[0] = Thi\n", "u[N] = Tlo\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0, 10.0 ]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_crank_nicolson_solve(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "plt.title(\"Heat equation\")\n", "plt.xlabel('${x}$ [cm]')\n", "plt.ylabel('${T}$ [Celsius]')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]) + \" s\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Heat equation in two dimensions\n", "\n", "In two dimensions the heat equation reads\n", "$$\n", "\\frac{\\partial u}{\\partial t} = D \\, \\left[ \\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2} \\right].\n", "$$\n", "\n", "This equation describes the time evolution of $u(t,x,y)$ given initial profile\n", "$$\n", "u(t=0,x,y) = u_0(x,y),\n", "$$\n", "and boundary conditions\n", "\\begin{align*}\n", "u(t,x=0,y) & = u_{\\rm left}(t;y), \\\\\n", "u(t,x=L,y) & = u_{\\rm right}(t;y), \\\\\n", "u(t,x=0,y) & = u_{\\rm bottom}(t;x), \\\\\n", "u(t,x=L,y) & = u_{\\rm top}(t;x).\n", "\\end{align*}\n", "\n", "Now we have two perform discretization in both $x$ and $y$ directions.\n", "Taking the same step size $a$ in both directions, we obtain the following discretized FTCS scheme:\n", "\n", "$$\n", "u^{n+1}_{i,j} = u^n_{i,j} + r \\, (u^n_{i+1,j} - 2u^n_{i,j} + u^n_{i-1,j})\n", "+ r \\, (u^n_{i,j+1} - 2u^n_{i,j} + u^n_{i,j-1}), \\qquad i = 1 \\ldots N-1, \\quad j = 1 \\ldots M-1.\n", "$$\n", "Here, as before,\n", "$$\n", "r \\equiv \\frac{Dh}{a^2},\n", "$$\n", "$N = L_x/a$, $M = L_y/a$, and\n", "$$\n", "u^n_{i,j} = u(t + hn, ai, aj).\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the 2D FTCS scheme in the time direction\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def heat_FTCS_iteration_2D(u, r):\n", " N, M = u.shape\n", " \n", " unew = np.empty_like(u)\n", " \n", " # Boundary conditions\n", " unew[ 0, :] = u[ 0, :]\n", " unew[N-1, :] = u[N-1, :]\n", " unew[ :, 0] = u[ :, 0]\n", " unew[ :, M-1] = u[ :, M-1]\n", " \n", " # FTCS scheme\n", " for i in range(1, M-1):\n", " for j in range(1, N-1):\n", " unew[i, j] = u[i, j] + r * (u[i+1, j] - 2 * u[i, j] + u[i-1, j]) + r * (u[i, j+1] - 2 * u[i, j] + u[i, j-1])\n", " \n", " return unew\n", "\n", "# Perform nsteps 2D FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def heat_FTCS_solve_2D(u0, h, nsteps, a, D = 1.):\n", " u = u0.copy()\n", " r = h * D / a**2\n", " for i in range(nsteps):\n", " u = heat_FTCS_iteration_2D(u, r)\n", " \n", " return u\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "M = N\n", "a = L/N # Grid spacing\n", "h = 1e-4 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with FTCS scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1,N+1],float)\n", "# Initial temperature\n", "u[:,:] = Tmid\n", "# Boundary conditions\n", "u[0,:] = Tlo\n", "u[N,:] = Tlo\n", "u[:,0] = Tlo\n", "u[:,M] = Thi\n", "\n", "# For the output\n", "times = [ 0.01, 0.1, 0.4, 1.0]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "yk = [k*a for k in range(0,M+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " u = heat_FTCS_solve_2D(u, h, nsteps, a, D)\n", " profiles.append(u.copy())\n", " current_time = time\n", " \n", "# Plot the initial and final temperature profiles\n", "fig, ax = plt.subplots(1, len(times), figsize=(12, 5))\n", "for i in range(len(ax)):\n", " ax[i].imshow(profiles[i].T, vmax=Thi, vmin=Tlo, origin=\"lower\", extent=[0,L,0,L])\n", " ax[i].set_title(\"t = \" + str(times[i]) + \" s\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "from IPython.display import clear_output, display\n", "from matplotlib.animation import FuncAnimation\n", "\n", "# Constants\n", "L = 0.01 # Thickness of steel in meters\n", "D = 4.25e-6 # Thermal diffusivity\n", "N = 100 # Number of divisions in grid\n", "M = N\n", "a = L/N # Grid spacing\n", "h = 1e-4 # Time-step (in s)\n", "\n", "print(\"Solving the heat equation with FTCS scheme\")\n", "print(\"r = h*D/a^2 =\",h*D/a**2)\n", "\n", "Tlo = 0.0 # Low temperature in Celsius\n", "Tmid = 20.0 # Intermediate temperature in Celsius\n", "Thi = 50.0 # High temperature in Celsius\n", "\n", "# Initialize\n", "u = np.zeros([N+1,N+1],float)\n", "# Initial temperature\n", "u[:,:] = Tmid\n", "# Boundary conditions\n", "u[0,:] = Tlo\n", "u[N,:] = Tlo\n", "u[:,0] = Tlo\n", "u[:,M] = Thi\n", "\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Heat equation\", fontsize = 18)\n", "# plt.subplots_adjust(top=1.)\n", "# function that draws each frame of the animation\n", "current_time = 0\n", "def animate_heat_equation(i):\n", " ax.clear()\n", " global fps\n", " niters = round((0.1/fps)/h)\n", " global u, current_time\n", " u = heat_FTCS_solve_2D(u, h, niters, a, D)\n", " current_time += h * niters\n", " clear_output(wait=True)\n", " print(\"t =\", current_time, \"/\", tend)\n", " \n", " ax.imshow(u.T, vmax=Thi, vmin=Tlo, origin=\"lower\", extent=[0,L,0,L])\n", " ax.set_title(\"t = \" + \"{:.3f}\".format(current_time) + \" s\")\n", "\n", "fps = 30\n", "tend = 2.\n", "ani = FuncAnimation(fig, animate_heat_equation, frames=fps * round(tend * 10), interval=1000/fps, repeat=False)\n", "# plt.show()\n", "\n", "ani.save(\"heat_equation_2D.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"heat_equation_2D.gif\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![heat2d](heat_equation_2D.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wave equation\n", "\n", "Wave equation is an example of a second-order linear PDE describing the waves and standing wave fields. In one dimensions it reads\n", "$$\n", "\\frac{\\partial^2 \\phi}{\\partial t^2} = v^2 \\frac{\\partial^2 \\phi}{\\partial x^2}.\n", "$$\n", "\n", "Since it is a 2nd order PDE, it is supplemented by initial conditions for both $\\phi(t=0,x)$ and $\\phi_t'(t=0,x)$:\n", "\\begin{align*}\n", "\\phi(t=0,x) & = \\phi_0(x) ,\\\\\n", "\\phi_t'(t=0,x) & = \\phi'_0(x).\n", "\\end{align*}\n", "\n", "The boundary conditions can be of either Dirichlet\n", "\\begin{align*}\n", "\\phi(t,x=0) & = \\phi_{\\rm left}(t) ,\\\\\n", "\\phi(t,x=L) & = \\phi_{\\rm right}(t),\n", "\\end{align*}\n", "or Neumann \n", "\\begin{align*}\n", "\\phi'_x(t,x=0) & = \\phi'_{\\rm left}(t) ,\\\\\n", "\\phi'_x(t,x=L) & = \\phi'_{\\rm right}(t),\n", "\\end{align*}\n", "forms.\n", "\n", "We shall focus on the Dirichlet form.\n", "\n", "## Finite difference approach\n", "\n", "\n", "To deal with the second-order time derivative we denote\n", "$$\n", "\\psi(t,x) \\equiv \\frac{\\partial \\phi}{\\partial t}.\n", "$$\n", "This way we are dealing with a system of first-order (in $t$) PDEs\n", "\\begin{align*}\n", "\\frac{\\partial \\phi}{\\partial t} = \\psi(t,x), \\\\\n", "\\frac{\\partial \\psi}{\\partial t} = v^2 \\frac{\\partial^2 \\phi}{\\partial x^2}.\n", "\\end{align*}\n", "\n", "To apply the finite difference method we first approximate the derivative $\\partial^2 \\phi / \\partial x^2$ by the lowest order central difference, just like for the heat equation,\n", "$$\n", "\\frac{\\partial^2 \\phi(t,x)}{\\partial x^2} \\approx \\frac{\\phi(t,x+a) - 2\\phi(t,x) + \\phi(t,x-a)}{a^2}.\n", "$$\n", "\n", "To solve the PDEs numerically we apply the same procedure as for the heat equation, but for $\\phi(t,x)$ and $\\psi(t,x)$ simultaneously.\n", "Denoting $\\phi(t = nh,x = ka) = \\phi^n_k$ and $\\psi(t = nh,x = ka) = \\psi^n_k$ we get\n", "\n", "### FTCS scheme\n", "\\begin{align*}\n", "\\phi^{n+1}_k & = \\phi^{n}_k + h \\psi^n_k, \\\\\n", "\\psi^{n+1}_k & = \\psi^n_k + r \\, (\\phi^n_{k+1} - 2\\phi^n_k + \\phi^n_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "\\end{align*}\n", "\n", "### Implicit scheme\n", "\\begin{align*}\n", "\\phi^{n+1}_k & = \\phi^{n}_k + h \\psi^{n+1}_k, \\\\\n", "\\psi^{n+1}_k & = \\psi^n_k + r \\, (\\phi^{n+1}_{k+1} - 2\\phi^{n+1}_k + \\phi^{n+1}_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "\\end{align*}\n", "\n", "### Crank-Nicolson scheme\n", "\\begin{align*}\n", "\\phi^{n+1}_k & = \\phi^{n}_k + \\frac{h}{2} \\left[ \\psi^{n+1}_k + \\psi^{n}_k\\right], \\\\\n", "\\psi^{n+1}_k & = \\psi^n_k + \\frac{r}{2} \\, (\\phi^{n+1}_{k+1} - 2\\phi^{n+1}_k + \\phi^{n+1}_{k-1}) + \\frac{r}{2} \\, (\\phi^{n}_{k+1} - 2\\phi^{n}_k + \\phi^{n}_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the FTCS scheme in the time direction\n", "# h is the time step\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def wave_FTCS_iteration(phi, psi, h, r):\n", " N = len(phi) - 1\n", " \n", " phinew = np.empty_like(phi)\n", " psinew = np.empty_like(psi)\n", " \n", " # Boundary conditions (here static Dirichlet)\n", " phinew[0] = phi[0]\n", " phinew[N] = phi[N]\n", " psinew[0] = 0.\n", " psinew[N] = 0.\n", " \n", " \n", " # FTCS scheme\n", " for i in range(1,N):\n", " phinew[i] = phi[i] + h * psi[i]\n", " psinew[i] = psi[i] + r * (phi[i+1] - 2 * phi[i] + phi[i-1])\n", " \n", " return phinew, psinew\n", "\n", "\n", "# Perform nsteps FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def wave_FTCS_solve(phi0, psi0, h, nsteps, a, v = 1.):\n", " phi = phi0.copy()\n", " psi = psi0.copy()\n", " r = h * v**2 / a**2\n", " for i in range(nsteps):\n", " phi, psi = wave_FTCS_iteration(phi, psi, h, r)\n", " \n", " return phi, psi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "v = 0.1 # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step\n", "\n", "print(\"Solving the wave equation with FTCS scheme\")\n", "print(\"r = h*v^2/a^2 =\",h*v**2/a**2)\n", "\n", "\n", "\n", "# Initialize\n", "phi = np.array([np.sin(k*np.pi/N) for k in range(N+1)])\n", "psi = np.zeros([N+1],float)\n", "\n", "# For the output\n", "times = [ 1., 5., 10., 18.]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " phi, psi = wave_FTCS_solve(phi, psi, h, nsteps, a, v)\n", " profiles.append(phi.copy())\n", " current_time = time\n", " \n", "plt.title(\"Wave equation with FTCS scheme\")\n", "plt.xlabel('${x}$')\n", "plt.ylabel('Displacement')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "from matplotlib.animation import FuncAnimation\n", "\n", "# Constants\n", "L = 1 # Length\n", "v = 0.1 # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step\n", "\n", "# For plotting\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "# Initialize for the three schemes\n", "phi1 = np.array([np.sin(k*np.pi/N) for k in range(N+1)])\n", "psi1 = np.zeros([N+1],float)\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Wave equation\", fontsize = 18)\n", "current_time = 0\n", "def animate_wave_equation_1D_FTCS(i):\n", " ax.clear()\n", " global fps, tottime, tend\n", " nsteps = round((tend/tottime/fps)/h)\n", " global phi1, psi1, phi2, psi2, phi3, psi3, current_time\n", " phi1, psi1 = wave_FTCS_solve(phi1, psi1, h, nsteps, a, v)\n", " current_time += h * nsteps\n", " # print(current_time)\n", " \n", " plt.title(\"${t=}$\" + \"{:.3f}\".format(current_time))\n", " plt.xlabel('${x}$')\n", " plt.ylabel('Displacement')\n", " plt.xlim(0,L)\n", " plt.ylim(-1.1,1.1)\n", " plt.plot(xk,phi1,label=\"FTCS\")\n", " plt.legend(loc=\"upper right\")\n", " # plt.show()\n", "\n", "fps = 30\n", "tend = 20.\n", "tottime = 10.\n", "ani = FuncAnimation(fig, animate_wave_equation_1D_FTCS, frames=round(fps * tottime), interval=1000/fps, repeat=False)\n", "#plt.show()\n", "\n", "ani.save(\"wave_equation_1D_FTCS.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"wave_equation_1D_FTCS.gif\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implicit scheme\n", "\\begin{align*}\n", "\\phi^{n+1}_k & = \\phi^{n}_k + h \\psi^{n+1}_k, \\\\\n", "\\psi^{n+1}_k & = \\psi^n_k + r \\, (\\phi^{n+1}_{k+1} - 2\\phi^{n+1}_k + \\phi^{n+1}_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "\\end{align*}\n", "\n", "Substituting the first equation into the second one gets the tridiagonal system of linear equations for $\\psi^{n+1}_k$:\n", "$$\n", "-rh \\psi^{n+1}_{k+1} + (1+2rh) \\psi^{n+1}_k - rh \\psi^{n+1}_{k-1} = \\psi^n_k + r \\, (\\phi^{n}_{k+1} - 2\\phi^{n}_k + \\phi^{n}_{k-1}), \\quad k = 1\\ldots N-1\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the FTCS scheme in the time direction\n", "# h is the time step\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def wave_implicit_iteration(phi, psi, h, r):\n", " N = len(phi) - 1\n", " \n", " phinew = np.empty_like(phi)\n", " psinew = np.empty_like(psi)\n", " \n", " # Boundary conditions (here static Dirichlet)\n", " phinew[0] = phi[0]\n", " phinew[N] = phi[N]\n", " psinew[0] = 0.\n", " psinew[N] = 0.\n", " \n", " # Tridiagonal system for psi\n", " d = np.full(N-1, 1+2.*r*h)\n", " ud = np.full(N-1, -r*h)\n", " ld = np.full(N-1, -r*h)\n", " v = np.array(psi[1:N] + r*phi[2:] - 2*r*phi[1:N] + r*phi[:-2])\n", " v[0] += r * h * psi[0]\n", " v[N-2] += r * h * psi[N]\n", " \n", " psinew[1:N] = linsolve_tridiagonal(d,ld,ud,v)\n", " \n", " # Final step for phi\n", " for i in range(1,N):\n", " phinew[i] = phi[i] + h * psinew[i]\n", " \n", " return phinew, psinew\n", "\n", "\n", "# Perform nsteps FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def wave_implicit_solve(phi0, psi0, h, nsteps, a, v = 1.):\n", " phi = phi0.copy()\n", " psi = psi0.copy()\n", " r = h * v**2 / a**2\n", " for i in range(nsteps):\n", " phi, psi = wave_implicit_iteration(phi, psi, h, r)\n", " \n", " return phi, psi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "v = 0.1 # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step\n", "\n", "print(\"Solving the wave equation with FTCS scheme\")\n", "print(\"r = h*v^2/a^2 =\",h*v**2/a**2)\n", "\n", "\n", "\n", "# Initialize\n", "phi = np.array([np.sin(k*np.pi/N) for k in range(N+1)])\n", "psi = np.zeros([N+1],float)\n", "\n", "# For the output\n", "times = [ 1., 5., 10., 20., 190., 290.]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " phi, psi = wave_implicit_solve(phi, psi, h, nsteps, a, v)\n", " profiles.append(phi.copy())\n", " current_time = time\n", " \n", "plt.title(\"Wave equation with implicit scheme\")\n", "plt.xlabel('${x}$')\n", "plt.ylabel('Displacement')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Crank-Nicolson scheme\n", "\\begin{align*}\n", "\\phi^{n+1}_k & = \\phi^{n}_k + \\frac{h}{2} \\left[ \\psi^{n+1}_k + \\psi^{n}_k\\right], \\\\\n", "\\psi^{n+1}_k & = \\psi^n_k + \\frac{r}{2} \\, (\\phi^{n+1}_{k+1} - 2\\phi^{n+1}_k + \\phi^{n+1}_{k-1}) + \\frac{r}{2} \\, (\\phi^{n}_{k+1} - 2\\phi^{n}_k + \\phi^{n}_{k-1}), \\qquad k = 1 \\ldots N-1.\n", "\\end{align*}\n", "\n", "Substituting the first equation into the second one gets the tridiagonal system of linear equations for $\\psi^{n+1}_k$:\n", "$$\n", "-rh \\psi^{n+1}_{k+1} + 2(1+rh) \\psi^{n+1}_k - rh \\psi^{n+1}_{k-1} = 2 \\psi^n_k + 2r \\, (\\phi^{n}_{k+1} - 2\\phi^{n}_k + \\phi^{n}_{k-1}) + rh \\, (\\psi^{n}_{k+1} - 2\\psi^{n}_k + \\psi^{n}_{k-1}), \\quad k = 1\\ldots N-1.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Single iteration of the FTCS scheme in the time direction\n", "# h is the time step\n", "# r = Dh/a^2 is the dimensionless parameter\n", "def wave_crank_nicolson_iteration(phi, psi, h, r):\n", " N = len(phi) - 1\n", " \n", " phinew = np.empty_like(phi)\n", " psinew = np.empty_like(psi)\n", " \n", " # Boundary conditions (here static Dirichlet)\n", " phinew[0] = phi[0]\n", " phinew[N] = phi[N]\n", " psinew[0] = 0.\n", " psinew[N] = 0.\n", " \n", " # Tridiagonal system for psi\n", " d = np.full(N-1, 2*(1+r*h))\n", " ud = np.full(N-1, -r*h)\n", " ld = np.full(N-1, -r*h)\n", " v = np.array(2*psi[1:N] + 2*r*phi[2:] - 4*r*phi[1:N] + 2*r*phi[:-2] \\\n", " + r*h*psi[2:] - 2*r*h*psi[1:N] + r*h*psi[:-2])\n", " v[0] += r * h * psi[0]\n", " v[N-2] += r * h * psi[N]\n", " \n", " psinew[1:N] = linsolve_tridiagonal(d,ld,ud,v)\n", " \n", " # Final step for phi\n", " for i in range(1,N):\n", " phinew[i] = phi[i] + h * (psinew[i] + psi[i]) / 2\n", " \n", " return phinew, psinew\n", "\n", "\n", "# Perform nsteps FTCS time iterations for the heat equation\n", "# u0: the initial profile\n", "# h: the size of the time step\n", "# nsteps: number of time steps\n", "# a: the spatial cell size\n", "# D: the diffusion constant\n", "def wave_crank_nicolson_solve(phi0, psi0, h, nsteps, a, v = 1.):\n", " phi = phi0.copy()\n", " psi = psi0.copy()\n", " r = h * v**2 / a**2\n", " for i in range(nsteps):\n", " phi, psi = wave_crank_nicolson_iteration(phi, psi, h, r)\n", " \n", " return phi, psi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "v = 0.1 # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step\n", "\n", "print(\"Solving the wave equation with FTCS scheme\")\n", "print(\"r = h*v^2/a^2 =\",h*v**2/a**2)\n", "\n", "\n", "\n", "# Initialize\n", "phi = np.array([np.sin(k*np.pi/N) for k in range(N+1)])\n", "psi = np.zeros([N+1],float)\n", "\n", "# For the output\n", "times = [ 1., 5., 10., 20., 190., 290.]\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " print(nsteps)\n", " phi, psi = wave_crank_nicolson_solve(phi, psi, h, nsteps, a, v)\n", " profiles.append(phi.copy())\n", " current_time = time\n", " \n", "plt.title(\"Wave equation with Crank-Nicolson scheme\")\n", "plt.xlabel('${x}$')\n", "plt.ylabel('Displacement')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the three" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "from matplotlib.animation import FuncAnimation\n", "\n", "# Constants\n", "L = 1 # Length\n", "v = 0.1 # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-2 # Time-step\n", "\n", "# For plotting\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "# Initialize for the three schemes\n", "phi1 = np.array([np.sin(k*np.pi/N) for k in range(N+1)])\n", "psi1 = np.zeros([N+1],float)\n", "phi2 = phi1.copy()\n", "psi2 = psi1.copy()\n", "phi3 = phi1.copy()\n", "psi3 = psi1.copy()\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Wave equation\", fontsize = 18)\n", "current_time = 0\n", "def animate_wave_equation_1D(i):\n", " ax.clear()\n", " global fps, tottime, tend\n", " nsteps = round((tend/tottime/fps)/h)\n", " global phi1, psi1, phi2, psi2, phi3, psi3, current_time\n", " phi1, psi1 = wave_FTCS_solve(phi1, psi1, h, nsteps, a, v)\n", " phi2, psi2 = wave_implicit_solve(phi2, psi2, h, nsteps, a, v)\n", " phi3, psi3 = wave_crank_nicolson_solve(phi3, psi3, h, nsteps, a, v)\n", " current_time += h * nsteps\n", " # print(current_time)\n", " \n", " plt.title(\"${t=}$\" + \"{:.3f}\".format(current_time))\n", " plt.xlabel('${x}$')\n", " plt.ylabel('Displacement')\n", " plt.xlim(0,L)\n", " plt.ylim(-1.1,1.1)\n", " plt.plot(xk,phi1,label=\"FTCS\")\n", " plt.plot(xk,phi2,label=\"Implicit\")\n", " plt.plot(xk,phi3,label=\"Crank-Nicolson\")\n", " plt.legend(loc=\"upper right\")\n", " # plt.show()\n", "\n", "fps = 30\n", "tend = 200.\n", "tottime = 10.\n", "ani = FuncAnimation(fig, animate_wave_equation_1D, frames=round(fps * tottime), interval=1000/fps, repeat=False)\n", "#plt.show()\n", "\n", "ani.save(\"wave_equation_1D.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"wave_equation_1D.gif\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exercise 9.5:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "d = 0.1 # Initial profile\n", "sig = 0.3 # Initial profile\n", "v = 100. # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-6 # Time-step\n", "\n", "print(\"Solving the wave equation with FTCS scheme\")\n", "print(\"r = h*v^2/a^2 =\",h*v**2/a**2)\n", "\n", "\n", "# Initialize\n", "phi = np.zeros([N+1],float)\n", "psi = np.array([a*k * (L - a*k) / L**2 * np.exp(-(k*a-d)**2/2/sig**2) for k in range(N+1)])\n", "\n", "# For the output\n", "times = [ 1.e-4, 2.e-4, 5.e-4, 8.e-4, 1.5e-3]\n", "for i in range(len(times)):\n", " times[i] *= 50.\n", "profiles = []\n", "xk = [k*a for k in range(0,N+1)]\n", "\n", "current_time = 0.\n", "for time in times:\n", " nsteps = round((time - current_time)/h)\n", " phi, psi = wave_FTCS_solve(phi, psi, h, nsteps, a, v)\n", " profiles.append(phi.copy())\n", " current_time = time\n", " \n", "plt.title(\"Wave equation with FTCS scheme\")\n", "plt.xlabel('${x}$')\n", "plt.ylabel('Displacement')\n", "for i in range(len(times)):\n", " plt.plot(xk,profiles[i],label=\"${t=}$\" + str(times[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "d = 0.1 # Initial profile\n", "sig = 0.3 # Initial profile\n", "v = 100. # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-6 # Time-step\n", "\n", "\n", "# Initialize\n", "phi1 = np.zeros([N+1],float)\n", "psi1 = np.array([a*k * (L - a*k) / L**2 * np.exp(-(k*a-d)**2/2/sig**2) for k in range(N+1)])\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Wave equation\", fontsize = 18)\n", "current_time = 0\n", "def animate_wave_equation_1D_pulse_FTCS(i):\n", " ax.clear()\n", " global fps, tottime, tend\n", " nsteps = round((tend/tottime/fps)/h)\n", " global phi1, psi1, current_time\n", " phi1, psi1 = wave_FTCS_solve(phi1, psi1, h, nsteps, a, v)\n", " current_time += h * nsteps\n", " # print(current_time)\n", " \n", " plt.title(\"${t=}$\" + \"{:.3f}\".format(current_time))\n", " plt.xlabel('${x}$')\n", " plt.ylabel('Displacement')\n", " plt.xlim(0,L)\n", " plt.ylim(-0.0005,0.0005)\n", " plt.plot(xk,phi1,label=\"FTCS\")\n", " plt.legend(loc=\"upper right\")\n", " # plt.show()\n", "\n", "fps = 30\n", "tend = 0.1\n", "tottime = 10.\n", "ani = FuncAnimation(fig, animate_wave_equation_1D_pulse_FTCS, frames=round(fps * tottime), interval=1000/fps, repeat=False)\n", "#plt.show()\n", "\n", "ani.save(\"wave_equation_1D_pulse_FTCS.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"wave_equation_1D_pulse_FTCS.gif\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare with implicit and Crank-Nicolson" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "# Constants\n", "L = 1 # Length\n", "d = 0.1 # Initial profile\n", "sig = 0.3 # Initial profile\n", "v = 100. # Wave propagation speed\n", "N = 100 # Number of divisions in grid\n", "a = L/N # Grid spacing\n", "h = 1e-6 # Time-step\n", "\n", "\n", "# Initialize\n", "phi1 = np.zeros([N+1],float)\n", "psi1 = np.array([a*k * (L - a*k) / L**2 * np.exp(-(k*a-d)**2/2/sig**2) for k in range(N+1)])\n", "phi2 = phi1.copy()\n", "psi2 = psi1.copy()\n", "phi3 = phi1.copy()\n", "psi3 = psi1.copy()\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "fig.set_size_inches(7, 5, forward=True)\n", "fig.suptitle(\"Wave equation\", fontsize = 18)\n", "current_time = 0\n", "def animate_wave_equation_1D_pulse(i):\n", " ax.clear()\n", " global fps, tottime, tend\n", " nsteps = round((tend/tottime/fps)/h)\n", " global phi1, psi1, phi2, psi2, phi3, psi3, current_time\n", " phi1, psi1 = wave_FTCS_solve(phi1, psi1, h, nsteps, a, v)\n", " phi2, psi2 = wave_implicit_solve(phi2, psi2, h, nsteps, a, v)\n", " phi3, psi3 = wave_crank_nicolson_solve(phi3, psi3, h, nsteps, a, v)\n", " current_time += h * nsteps\n", " # print(current_time)\n", " \n", " plt.title(\"${t=}$\" + \"{:.3f}\".format(current_time))\n", " plt.xlabel('${x}$')\n", " plt.ylabel('Displacement')\n", " plt.xlim(0,L)\n", " plt.ylim(-0.0005,0.0005)\n", " plt.plot(xk,phi1,label=\"FTCS\")\n", " plt.plot(xk,phi2,label=\"Implicit\")\n", " plt.plot(xk,phi3,label=\"Crank-Nicolson\")\n", " plt.legend(loc=\"upper right\")\n", " # plt.show()\n", "\n", "fps = 30\n", "tend = 0.2\n", "tottime = 10.\n", "ani = FuncAnimation(fig, animate_wave_equation_1D_pulse, frames=round(fps * tottime), interval=1000/fps, repeat=False)\n", "#plt.show()\n", "\n", "ani.save(\"wave_equation_1D_pulse.gif\")\n", "\n", "plt.close(fig)\n", "from IPython.display import display, Image, clear_output\n", "display(Image(filename=\"wave_equation_1D_pulse.gif\"))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.13" } }, "nbformat": 4, "nbformat_minor": 4 }