{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Multi-dimensional nonlinear PDE problems\n", "\n", "The fundamental ideas in the\n", "derivation of $F_i$ and $J_{i,j}$ in the 1D model problem\n", "are easily generalized to multi-dimensional problems.\n", "Nevertheless, the expressions involved are slightly different, with\n", "derivatives in $x$ replaced by $\\nabla$, so we present some\n", "examples below in detail.\n", "\n", "\n", "## Finite difference discretization\n", "
\n", "\n", "A typical diffusion equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_t = \\nabla\\cdot(\\dfc(u)\\nabla u) + f(u),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "can be discretized by (e.g.) a Backward Euler scheme,\n", "which in 2D can be written" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t^- u = D_x\\overline{\\dfc(u)}^xD_x u\n", "+ D_y\\overline{\\dfc(u)}^yD_y u + f(u)]_{i,j}^n\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do not dive into the details of handling boundary conditions now.\n", "Dirichlet and Neumann conditions are handled as in\n", "corresponding linear, variable-coefficient diffusion problems.\n", "\n", "Writing the scheme out, putting the unknown values on the\n", "left-hand side and known values on the right-hand side, and\n", "introducing $\\Delta x=\\Delta y=h$ to save some writing, one gets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u^n_{i,j} &- \\frac{\\Delta t}{h^2}(\n", " \\frac{1}{2}(\\dfc(u_{i,j}^n) + \\dfc(u_{i+1,j}^n))(u_{i+1,j}^n-u_{i,j}^n)\\\\ \n", "&\\quad -\n", "\\frac{1}{2}(\\dfc(u_{i-1,j}^n) + \\dfc(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j}^n) \\\\ \n", "&\\quad +\n", " \\frac{1}{2}(\\dfc(u_{i,j}^n) + \\dfc(u_{i,j+1}^n))(u_{i,j+1}^n-u_{i,j}^n)\\\\ \n", "&\\quad -\n", " \\frac{1}{2}(\\dfc(u_{i,j-1}^n) + \\dfc(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j-1}^n))\n", "- \\Delta tf(u_{i,j}^n) = u^{n-1}_{i,j}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This defines a nonlinear algebraic system on the form $A(u)u=b(u)$.\n", "\n", "### Picard iteration\n", "\n", "The most recently computed values $u^{-}$ of $u^n$ can be\n", "used in $\\dfc$ and $f$ for a Picard iteration,\n", "or equivalently, we solve $A(u^{-})u=b(u^{-})$.\n", "The result is a linear system of the same type as arising\n", "from $u_t = \\nabla\\cdot (\\dfc(\\x)\\nabla u) + f(\\x,t)$.\n", "\n", "The Picard iteration scheme can also be expressed in operator notation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t^- u = D_x\\overline{\\dfc(u^{-})}^xD_x u\n", "+ D_y\\overline{\\dfc(u^{-})}^yD_y u + f(u^{-})]_{i,j}^n\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Newton's method\n", "\n", "As always, Newton's method is technically more involved than\n", "Picard iteration. We first define\n", "the nonlinear algebraic equations to be solved, drop the\n", "superscript $n$ (use $u$ for $u^n$), and introduce $u^{(1)}$ for $u^{n-1}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "F_{i,j} &= u_{i,j} - \\frac{\\Delta t}{h^2}(\\\\ \n", "&\\qquad \\frac{1}{2}(\\dfc(u_{i,j}) + \\dfc(u_{i+1,j}))(u_{i+1,j}-u_{i,j}) -\\\\ \n", "&\\qquad \\frac{1}{2}(\\dfc(u_{i-1,j}) + \\dfc(u_{i,j}))(u_{i,j}-u_{i-1,j}) + \\\\ \n", "&\\qquad \\frac{1}{2}(\\dfc(u_{i,j}) + \\dfc(u_{i,j+1}))(u_{i,j+1}-u_{i,j}) -\\\\ \n", "&\\qquad \\frac{1}{2}(\\dfc(u_{i,j-1}) + \\dfc(u_{i,j}))(u_{i,j}-u_{i-1,j-1})) -\n", "\\Delta t\\, f(u_{i,j}) - u^{(1)}_{i,j} = 0\\thinspace .\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mathcal{I}_t is convenient to work with two indices $i$ and $j$ in 2D\n", "finite difference discretizations, but it complicates\n", "the derivation of the Jacobian, which then gets four indices.\n", "(Make sure you really understand the 1D version of this problem\n", "as treated in the section [nonlin:alglevel:1D:fd](#nonlin:alglevel:1D:fd).)\n", "The left-hand expression of an equation $F_{i,j}=0$ is to be\n", "differentiated with respect to each of the unknowns $u_{r,s}$\n", "(recall that this is short notation for $u_{r,s}^n$), $r\\in\\mathcal{I}_x$, $s\\in\\Iy$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "J_{i,j,r,s} = \\frac{\\partial F_{i,j}}{\\partial u_{r,s}}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Newton system to be solved in each iteration can be written as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sum_{r\\in\\mathcal{I}_x}\\sum_{s\\in\\Iy}J_{i,j,r,s}\\delta u_{r,s} = -F_{i,j},\n", "\\quad i\\in\\mathcal{I}_x,\\ j\\in\\Iy\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given $i$ and $j$, only a few $r$ and $s$ indices give nonzero\n", "contribution to the Jacobian since $F_{i,j}$ contains $u_{i\\pm 1,j}$,\n", "$u_{i,j\\pm 1}$, and $u_{i,j}$. This means that $J_{i,j,r,s}$ has\n", "nonzero contributions only if $r=i\\pm 1$, $s=j\\pm 1$, as well as $r=i$\n", "and $s=j$. The corresponding terms in $J_{i,j,r,s}$ are\n", "$J_{i,j,i-1,j}$, $J_{i,j,i+1,j}$, $J_{i,j,i,j-1}$, $J_{i,j,i,j+1}$\n", "and $J_{i,j,i,j}$. Therefore, the left-hand side of the Newton\n", "system, $\\sum_r\\sum_s J_{i,j,r,s}\\delta u_{r,s}$ collapses to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", " J_{i,j,r,s}\\delta u_{r,s} =\n", "J_{i,j,i,j}\\delta u_{i,j} & +\n", "J_{i,j,i-1,j}\\delta u_{i-1,j} +\n", "J_{i,j,i+1,j}\\delta u_{i+1,j} +\n", "J_{i,j,i,j-1}\\delta u_{i,j-1}\\\\ \n", "& +\n", "J_{i,j,i,j+1}\\delta u_{i,j+1}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The specific derivatives become" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "J_{i,j,i-1,j} &= \\frac{\\partial F_{i,j}}{\\partial u_{i-1,j}}\\\\ \n", "&= \\frac{\\Delta t}{h^2}(\\dfc^{\\prime}(u_{i-1,j})(u_{i,j}-u_{i-1,j})\n", "+ \\dfc(u_{i-1,j})(-1)),\\\\ \n", "J_{i,j,i+1,j} &= \\frac{\\partial F_{i,j}}{\\partial u_{i+1,j}}\\\\ \n", "&= \\frac{\\Delta t}{h^2}(-\\dfc^{\\prime}(u_{i+1,j})(u_{i+1,j}-u_{i,j})\n", "- \\dfc(u_{i-1,j})),\\\\ \n", "J_{i,j,i,j-1} &= \\frac{\\partial F_{i,j}}{\\partial u_{i,j-1}}\\\\ \n", "&= \\frac{\\Delta t}{h^2}(\\dfc^{\\prime}(u_{i,j-1})(u_{i,j}-u_{i,j-1})\n", "+ \\dfc(u_{i,j-1})(-1)),\\\\ \n", "J_{i,j,i,j+1} &= \\frac{\\partial F_{i,j}}{\\partial u_{i,j+1}}\\\\ \n", "&= \\frac{\\Delta t}{h^2}(-\\dfc^{\\prime}(u_{i,j+1})(u_{i,j+1}-u_{i,j})\n", "- \\dfc(u_{i,j-1}))\\thinspace .\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $J_{i,j,i,j}$ entry has a few more terms and is left as an\n", "exercise.\n", "Inserting the most recent approximation\n", "$u^{-}$ for $u$ in the $J$ and $F$ formulas and then\n", "forming $J\\delta u=-F$ gives the linear system to be solved\n", "in each Newton iteration. Boundary conditions will affect the\n", "formulas when any of the indices coincide with a boundary value\n", "of an index.\n", "\n", "\n", "## Continuation methods\n", "\n", "\n", "Picard iteration or Newton's method may diverge when solving PDEs with\n", "severe nonlinearities. Relaxation with $\\omega <1$\n", "may help, but in highly nonlinear problems it can be\n", "necessary to introduce a *continuation parameter* $\\Lambda$ in\n", "the problem: $\\Lambda =0$ gives a version of the\n", "problem that is easy to solve, while\n", "$\\Lambda =1$ is the target problem. The idea is then\n", "to increase $\\Lambda$ in steps, $\\Lambda_0=0 ,\\Lambda_1 <\\cdots <\\Lambda_n=1$,\n", "and use the solution from the problem with $\\Lambda_{i-1}$ as\n", "initial guess for the iterations in the problem corresponding\n", "to $\\Lambda_i$.\n", "\n", "The continuation method is easiest to understand through an example.\n", "Suppose we intend to solve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-\\nabla\\cdot\\left( ||\\nabla u||^q\\nabla u\\right) = f,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is an equation modeling the flow of a non-Newtonian fluid through\n", "a channel or pipe. For $q=0$ we have the Poisson equation (corresponding\n", "to a Newtonian fluid) and the problem is linear. A typical\n", "value for pseudo-plastic fluids may be $q_n=-0.8$. We can introduce\n", "the continuation parameter $\\Lambda\\in [0,1]$ such that\n", "$q=q_n\\Lambda$. Let $\\{\\Lambda_\\ell\\}_{\\ell=0}^n$ be the sequence of\n", "$\\Lambda$ values in $[0,1]$, with corresponding $q$ values\n", "$\\{q_\\ell\\}_{\\ell=0}^n$. We can then solve a sequence of problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-\\nabla\\cdot\\left( ||\\nabla u^\\ell||^q_\\ell\\nabla u^\\ell\\right) = f,\\quad\n", "\\ell = 0,\\ldots,n,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the initial guess for iterating on $u^{\\ell}$ is the\n", "previously computed solution $u^{\\ell-1}$. If a particular $\\Lambda_\\ell$\n", "leads to convergence problems, one may try a smaller\n", "increase in $\\Lambda$:\n", "$\\Lambda_* = \\frac{1}{2} (\\Lambda_{\\ell-1}+\\Lambda_\\ell)$,\n", "and repeat halving the step in $\\Lambda$ until convergence is reestablished." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 4 }