draft of complete parabolic module

content/pdes/parabolic.ipynb

@@ -35,6 +35,11 @@
     "\n",
     "To solve the heat equation for a one-dimensional domain over $0 \\leq x \\leq L$, we will need both initial conditions at $t = 0$ and boundary conditions at $x=0$ and $x=L$ (for all time). In terms of our nodal values, this means we need $T_i^{k=1}$ for $i = 1 \\ldots n$, where $n$ is the number of points, as well as information about $T_1^k$ and $T_n^k$ for all times $k$.\n",
     "\n",
+    "Based on the choices of finite differences, we have some information on the accuracy of this approach:\n",
+    "\n",
+    "- It is **second-order accurate** in space, with respect to $\\Delta x$\n",
+    "- It is **first-order accurate** in time, with respect to $\\Delta t$.\n",
+    "\n",
     "We can rearrange the above equation to obtain our recursion formula:\n",
     "\\begin{equation}\n",
     "T_i^{k+1} = \\left( T_{i+1}^k + T_{i-1}^k \\right) \\frac{\\alpha \\Delta t}{\\Delta x^2} + T_i^k  \\left( 1 - 2 \\frac{\\alpha \\Delta t}{\\Delta x^2} \\right) \\;.\n",
@@ -67,9 +72,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time step size: 0.011\n"
+     ]
+    }
+   ],
    "source": [
     "clear all\n",
     "\n",
@@ -81,6 +94,8 @@
     "% then choose the time step based on the Fourier number\n",
     "dt = Fo * dx^2 / alpha;\n",
     "\n",
+    "fprintf('Time step size: %5.3f\\n', dt);\n",
+    "\n",
     "x = [0 : dx : 1]; n = length(x);\n",
     "t = [0 : dt : 1]; m = length(t);\n",
     "\n",
@@ -152,9 +167,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 3,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time step size: 0.033\n"
+     ]
+    }
+   ],
    "source": [
     "clear all\n",
     "\n",
@@ -165,6 +188,8 @@
     "Fo = 0.75;\n",
     "dt = Fo * dx^2 / alpha;\n",
     "\n",
+    "fprintf('Time step size: %5.3f\\n', dt);\n",
+    "\n",
     "x = [0 : dx : 1]; n = length(x);\n",
     "t = [0 : dt : 1]; m = length(t);\n",
     "\n",
@@ -205,7 +230,7 @@
     "    [imind,cm] = rgb2ind(im,256); \n",
     "    % Write to GIF\n",
     "    if i == 1 \n",
-    "        imwrite(imind,cm,filename,'gif', 'Loopcount',2, 'DelayTime',1e-3); \n",
+    "        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3); \n",
     "    else \n",
     "        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3); \n",
     "    end\n",
@@ -230,7 +255,7 @@
    "source": [
     "In this case, the solution becomes unstable and blows up, leading to unphysical results.\n",
     "\n",
-    "For this explicit scheme, the choice of $\\Delta t$ is limited by the stability criterion. This means that we may be stuck using a small time-step size.\n",
+    "For this explicit scheme, which is **conditionally stable**, the choice of $\\Delta t$ is limited by the stability criterion. This means that we may be stuck using a small time-step size.\n",
     "\n",
     "Rather than being forced to use a very small time-step size, we can also explore *implicit schemes* that are unconditionally stable."
    ]
@@ -240,14 +265,47 @@
    "metadata": {},
    "source": [
     "## Implicit scheme\n",
-    "\n"
+    "\n",
+    "To solve the problem without being constrained by the stability criterion (and the associated restriction on time-step size), we can develop an *implicit* scheme for solving the 1D unsteady heat equation by evaluating the second derivative at time $t^{k+1}$, rather than at time $t^k$ as we did above:\n",
+    "\\begin{equation}\n",
+    "\\frac{T_i^{k+1} - T_i^k}{\\Delta t} = \\alpha \\frac{T_{i-1}^{k+1} - 2 T_i^{k+1} + T_{i+1}^{k+1}}{\\Delta x^2}\n",
+    "\\end{equation}\n",
+    "\n",
+    "Then, rearrange this to obtain a recursion formula with unknowns on the left-hand side and knowns on the right-hand side:\n",
+    "\\begin{equation}\n",
+    "T_{i-1}^{k+1} (\\text{Fo}) + T_i^{k+1} \\left(-2 \\text{Fo} - 1 \\right) + T_{i+1}^{k+1} (\\text{Fo}) = -T_i^k\n",
+    "\\end{equation}\n",
+    "where $\\text{Fo} = \\alpha \\Delta t / \\Delta x^2$ is the Fourier number.\n",
+    "\n",
+    "<figure>\n",
+    "  <center>\n",
+    "  <img src=\"../images/parabolic-implicit-stencil.png\" alt=\"stencil for implicit parabolic solution\" style=\"width: 350px;\"/>\n",
+    "  <figcaption>Figure: Stencil for implicit solution to heat equation</figcaption>\n",
+    "  </center>\n",
+    "</figure>\n",
+    "\n",
+    "Unlike the explicit scheme, with the implicit scheme we do not have a simple recursion formula that can be applied to calculate the next temperature value at each location. Instead, we need to solve a system of linear equations at each time step, to simultaneously get all the values of temperature at the next time step.\n",
+    "In other words, at each time step, we need to solve\n",
+    "\\begin{equation}\n",
+    "A \\mathbf{T}^{k+1} = - \\mathbf{T}^k = \\mathbf{b} \\;.\n",
+    "\\end{equation}\n",
+    "\n",
+    "Let's implement this to solve the same example:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 4,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time step size: 0.033\n"
+     ]
+    }
+   ],
    "source": [
     "clear all\n",
     "\n",
@@ -259,6 +317,8 @@
     "Fo = 0.75;\n",
     "dt = Fo * dx^2 / alpha;\n",
     "\n",
+    "fprintf('Time step size: %5.3f\\n', dt);\n",
+    "\n",
     "t = [0 : dt : 1]; m = length(t);\n",
     "\n",
     "T = zeros(m, n);\n",
@@ -297,6 +357,10 @@
     "end\n",
     "close\n",
     "\n",
+    "%% If you are working interactively, you can use this to make a movie in Matlab\n",
+    "%fig = figure;\n",
+    "%movie(fig, F, 2)\n",
+    "\n",
     "%% This generates a GIF of the results (for use in Jupyter Notebook)\n",
     "filename = 'parabolic_implicit_animated.gif';\n",
     "for i = 1 : length(F)\n",
@@ -323,6 +387,172 @@
     "</figure>"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With this implicit scheme, we can now use larger time-step sizes (resulting in values of the Fourier number greater than 0.5), because the method is **unconditionally stable**. Like the explicit scheme above, this approach is:\n",
+    "\n",
+    "- Second-order accurate in space, with respect to $\\Delta x$\n",
+    "- First-order accurate in time, with respect to $\\Delta t$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Crank-Nicolson scheme\n",
+    "\n",
+    "So far we have two methods to solve parabolic PDEs:\n",
+    "\n",
+    "- The explicit scheme, which is conditionally stable and first-order accurate in time\n",
+    "- The implicit scheme, which is unconditionally stable but also first-order accurate in time\n",
+    "\n",
+    "We might also want a scheme that is more accurate with respect to the time-step size $\\Delta t$. The **Crank-Nicolson method** is 2nd-order accurate in time and also unconditionally stable (because it is implicit). This method is essentially an average of the first-order accurate explicit and implicit methods:\n",
+    "\\begin{equation}\n",
+    "\\frac{T_i^{k+1} - T_i^k}{\\Delta t} = \\frac{\\alpha}{2} \\frac{T_{i-1}^{k} - 2 T_i^{k} + T_{i+1}^{k}}{\\Delta x^2} + \\frac{\\alpha}{2} \\frac{T_{i-1}^{k+1} - 2 T_i^{k+1} + T_{i+1}^{k+1}}{\\Delta x^2}\n",
+    "\\end{equation}\n",
+    "This method is based on using the average of the time derivative at the current and next time steps.\n",
+    "\n",
+    "<figure>\n",
+    "  <center>\n",
+    "  <img src=\"../images/parabolic-CN-stencil.png\" alt=\"stencil for Crank-Nicolson parabolic solution\" style=\"width: 350px;\"/>\n",
+    "  <figcaption>Figure: Stencil for Crank-Nicolson solution to heat equation</figcaption>\n",
+    "  </center>\n",
+    "</figure>\n",
+    "\n",
+    "If we rearrange that equation, we can get our recursion formula:\n",
+    "\\begin{equation}\n",
+    "T_{i-1}^{k+1} \\left( \\frac{\\text{Fo}}{2} \\right) + T_{i}^{k+1} \\left( Fo - 1 \\right) + T_{i+1}^{k+1} \\left( \\frac{\\text{Fo}}{2} \\right) = T_{i-1}^{k} \\left( -\\frac{\\text{Fo}}{2} \\right) + T_{i}^{k} \\left( \\text{Fo} - 1 \\right) + T_{i-1}^{k} \\left( -\\frac{\\text{Fo}}{2} \\right) \\;,\n",
+    "\\end{equation}\n",
+    "where $\\text{Fo} = \\alpha \\Delta t / \\Delta x^2$.\n",
+    "Like the first-order implicit scheme, advancing this solution in time requires solving a system of linear equations at each time step:\n",
+    "\\begin{equation}\n",
+    "A \\mathbf{T}^{k+1} = B \\mathbf{T}^k = \\mathbf{b} \\;,\n",
+    "\\end{equation}\n",
+    "where the coefficient matrices $A$ and $B$ will look like\n",
+    "\\begin{equation}\n",
+    "A = \\begin{bmatrix}\n",
+    "1 \\\\\n",
+    "\\frac{\\text{Fo}}{2} & (-\\text{Fo}-1) & \\frac{\\text{Fo}}{2} \\\\\n",
+    " & \\frac{\\text{Fo}}{2} & (-\\text{Fo}-1) & \\frac{\\text{Fo}}{2} \\\\\n",
+    " & & \\ddots & \\ddots & \\ddots & \\\\\n",
+    " & &        &        &        & 1\n",
+    "\\end{bmatrix}\n",
+    "\\quad\n",
+    "B = \\begin{bmatrix}\n",
+    "1 \\\\\n",
+    "\\frac{-\\text{Fo}}{2} & (\\text{Fo}-1) & \\frac{\\text{-Fo}}{2} \\\\\n",
+    " & \\frac{-\\text{Fo}}{2} & (\\text{Fo}-1) & \\frac{-\\text{Fo}}{2} \\\\\n",
+    " & & \\ddots & \\ddots & \\ddots & \\\\\n",
+    " & &        &        &        & 1\n",
+    "\\end{bmatrix}\n",
+    "\\end{equation}\n",
+    "The first and last rows of these matrices may differ, depending on the boundary conditions.\n",
+    "\n",
+    "Let's implement this method:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time step size: 0.033\n"
+     ]
+    }
+   ],
+   "source": [
+    "clear all\n",
+    "\n",
+    "alpha = 2.3e-1;\n",
+    "dx = 0.1;\n",
+    "x = [0 : dx : 1]; n = length(x);\n",
+    "\n",
+    "% choose a Fourier number that is deliberately past the explicit method stability limit\n",
+    "Fo = 0.75;\n",
+    "dt = Fo * dx^2 / alpha;\n",
+    "\n",
+    "fprintf('Time step size: %5.3f\\n', dt);\n",
+    "\n",
+    "t = [0 : dt : 1]; m = length(t);\n",
+    "\n",
+    "T = zeros(m, n);\n",
+    "\n",
+    "% Initial conditions\n",
+    "T(1,:) = 200;\n",
+    "\n",
+    "% boundary conditions\n",
+    "T(:,1) = 50;\n",
+    "T(:,n) = 50;\n",
+    "\n",
+    "plot(x, T(1,:))\n",
+    "axis([0 1 50 200]);\n",
+    "xlabel('Distance'); ylabel('Temperature');\n",
+    "F(1) = getframe(gcf);\n",
+    "\n",
+    "for k = 1 : m - 1\n",
+    "    A = zeros(n,n);\n",
+    "    B = zeros(n,n);\n",
+    "    for i = 1 : n\n",
+    "        if i == 1\n",
+    "            A(1,1) = 1;\n",
+    "            B(1,1) = 1;\n",
+    "        elseif i == n\n",
+    "            A(n,n) = 1;\n",
+    "            B(n,n) = 1;\n",
+    "        else\n",
+    "            A(i,i-1) = Fo/2;\n",
+    "            A(i,i) = -Fo - 1;\n",
+    "            A(i,i+1) = Fo/2;\n",
+    "            B(i,i-1) = -Fo/2;\n",
+    "            B(i,i) = Fo - 1;\n",
+    "            B(i,i+1) = -Fo/2;\n",
+    "        end\n",
+    "    end\n",
+    "    b = B * T(k,:)';\n",
+    "    T(k+1, :) = A \\ b;\n",
+    "    plot(x, T(k+1,:))\n",
+    "    axis([0 1 50 200]);\n",
+    "    xlabel('Distance'); ylabel('Temperature');\n",
+    "    F(k+1) = getframe(gcf);\n",
+    "end\n",
+    "close\n",
+    "\n",
+    "%% If you are working interactively, you can use this to make a movie in Matlab\n",
+    "%fig = figure;\n",
+    "%movie(fig, F, 2)\n",
+    "\n",
+    "%% This generates a GIF of the results (for use in Jupyter Notebook)\n",
+    "filename = 'parabolic_cranknicolson_animated.gif';\n",
+    "for i = 1 : length(F)\n",
+    "    im = frame2im(F(i)); \n",
+    "    [imind,cm] = rgb2ind(im,256); \n",
+    "    % Write to GIF\n",
+    "    if i == 1 \n",
+    "        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3); \n",
+    "    else \n",
+    "        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3); \n",
+    "    end\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<figure>\n",
+    "  <center>\n",
+    "  <img src=\"parabolic_cranknicolson_animated.gif\" alt=\"movie of Crank-Nicolson solution to parabolic PDE\" style=\"width: 500px;\"/>\n",
+    "  <figcaption>Figure: Solution to 1D heat equation with Crank-Nicolson method. Fo = 0.75</figcaption>\n",
+    "  </center>\n",
+    "</figure>"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,