8- nonlin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dealing with nonlinearity - ODEs\n",
    "<div id=\"nonlin:timediscrete:logistic\"></div>\n",
    "\n",
    "## Linear versus nonlinear equations\n",
    "\n",
    "### Algebraic equations\n",
    "\n",
    "A linear, scalar, algebraic equation in $x$ has the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "ax + b = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for arbitrary real constants $a$ and $b$. The unknown is a number $x$.\n",
    "All other algebraic equations, e.g., $x^2 + ax + b = 0$, are nonlinear.\n",
    "The typical feature in a nonlinear algebraic equation is that the unknown\n",
    "appears in products with itself, like $x^2$ or\n",
    "in functions that are infinite sums of products, like $e^x = 1 + x +\\frac{1}{2} x^2 +\n",
    "\\frac{1}{3!}x^3 + \\cdots$.\n",
    "\n",
    "We know how to solve a linear algebraic equation, $x=-b/a$, but there are\n",
    "no general closed formulas for finding the exact solutions of\n",
    "nonlinear algebraic equations, except for very special cases (quadratic\n",
    "equations constitute a primary example). A nonlinear algebraic equation\n",
    "may have no solution, one solution, or many solutions. The tools for\n",
    "solving nonlinear algebraic equations are *iterative methods*, where\n",
    "we construct a series of linear equations, which we know how to solve,\n",
    "and hope that the solutions of the linear equations converge to a\n",
    "solution of the nonlinear equation we want to solve.\n",
    "Typical methods for nonlinear algebraic equation equations are\n",
    "Newton's method, the Bisection method, and the Secant method.\n",
    "\n",
    "### Differential equations\n",
    "\n",
    "The unknown in a differential equation is a function and not a number.\n",
    "In a linear differential equation, all terms involving the unknown function\n",
    "are linear in the unknown function or its derivatives. Linear here means that\n",
    "the unknown function, or a derivative of it, is multiplied by a number or\n",
    "a known function. All other differential equations are non-linear.\n",
    "\n",
    "\n",
    "The easiest way to see if an equation is nonlinear, is to spot nonlinear terms\n",
    "where the unknown function or its derivatives are multiplied by\n",
    "each other. For example, in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{\\prime}(t) = -a(t)u(t) + b(t),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the terms involving the unknown function $u$ are linear: $u^{\\prime}$ contains\n",
    "the derivative of the unknown function multiplied by unity, and $au$ contains\n",
    "the unknown function multiplied by a known function.\n",
    "However,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{\\prime}(t) = u(t)(1 - u(t)),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is nonlinear because of the term $-u^2$ where the unknown function is\n",
    "multiplied by itself. Also"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is nonlinear because of the term $uu_x$ where the unknown\n",
    "function appears in a product with its derivative.\n",
    "(Note here that we use different notations for derivatives: $u^{\\prime}$\n",
    "or $du/dt$ for a function $u(t)$ of one variable,\n",
    "$\\frac{\\partial u}{\\partial t}$ or $u_t$ for a function of more than one\n",
    "variable.)\n",
    "\n",
    "Another example of a nonlinear equation is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{\\prime\\prime} + \\sin(u) =0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "because $\\sin(u)$ contains products of $u$, which becomes clear\n",
    "if we expand the function in a Taylor series:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\sin(u) = u - \\frac{1}{3} u^3 + \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Mathematical proof of linearity.**\n",
    "\n",
    "To really prove mathematically that some differential equation\n",
    "in an unknown $u$ is linear,\n",
    "show for each term $T(u)$ that with $u = au_1 + bu_2$ for\n",
    "constants $a$ and $b$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T(au_1 + bu_2) = aT(u_1) + bT(u_2){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the term $T(u) = (\\sin^2 t)u'(t)$ is linear because"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "T(au_1 + bu_2) &= (\\sin^2 t)(au_1(t) + b u_2(t))'\\\\\n",
    "& = a(\\sin^2 t)u_1'(t) + b(\\sin^2 t)u_2'(t)\\\\\n",
    "& =aT(u_1) + bT(u_2){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, $T(u)=\\sin u$ is nonlinear because"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T(au_1 + bu_2) = \\sin (au_1 + bu_2) \\neq a\\sin u_1 + b\\sin u_2{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple model problem\n",
    "\n",
    "A series of forthcoming examples will explain how to tackle\n",
    "nonlinear differential equations with various techniques.\n",
    "We start with the (scaled) logistic equation as model problem:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:eq\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u^{\\prime}(t) = u(t)(1 - u(t)) {\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:eq} \\tag{1}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a nonlinear ordinary differential equation (ODE)\n",
    "which will be solved by\n",
    "different strategies in the following.\n",
    "Depending on the chosen\n",
    "time discretization of ([1](#nonlin:timediscrete:logistic:eq)),\n",
    "the mathematical problem to be solved at every time level will\n",
    "either be a linear algebraic equation or a nonlinear\n",
    "algebraic equation.\n",
    "In the former case, the time discretization method transforms\n",
    "the nonlinear ODE into linear subproblems at each time level, and\n",
    "the solution is straightforward to find since linear algebraic equations\n",
    "are easy to solve. However,\n",
    "when the time discretization leads to nonlinear algebraic equations, we\n",
    "cannot (except in very rare cases) solve these without turning to\n",
    "approximate, iterative solution methods.\n",
    "\n",
    "The next subsections introduce various methods\n",
    "for solving nonlinear differential equations,\n",
    "using ([1](#nonlin:timediscrete:logistic:eq)) as model. We shall go through\n",
    "the following set of cases:\n",
    "\n",
    " * explicit time discretization methods (with no need to\n",
    "   solve nonlinear algebraic equations)\n",
    "\n",
    " * implicit Backward Euler time discretization, leading to nonlinear\n",
    "   algebraic equations solved by\n",
    "\n",
    "  * an exact analytical technique\n",
    "\n",
    "  * Picard iteration based on manual linearization\n",
    "\n",
    "  * a single Picard step\n",
    "\n",
    "  * Newton's method\n",
    "\n",
    "\n",
    " * implicit Crank-Nicolson time discretization and linearization\n",
    "   via a geometric mean formula\n",
    "\n",
    "Thereafter, we compare the performance of the various approaches. Despite\n",
    "the simplicity of ([1](#nonlin:timediscrete:logistic:eq)), the conclusions\n",
    "reveal typical features of the various methods in much more complicated\n",
    "nonlinear PDE problems.\n",
    "\n",
    "\n",
    "## Linearization by explicit time discretization\n",
    "<div id=\"nonlin:timediscrete:logistic:FE\"></div>\n",
    "\n",
    "\n",
    "Time discretization methods are divided into explicit and implicit\n",
    "methods. Explicit methods lead to a closed-form formula for\n",
    "finding new values of the unknowns, while implicit methods give\n",
    "a linear or nonlinear system of equations that couples (all) the\n",
    "unknowns at a new time level. Here we shall demonstrate that\n",
    "explicit methods may constitute an efficient way to deal with nonlinear\n",
    "differential equations.\n",
    "\n",
    "The Forward Euler\n",
    "method is an explicit method. When applied to\n",
    "([1](#nonlin:timediscrete:logistic:eq)), sampled at $t=t_n$, it results in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n+1} - u^n}{\\Delta t} = u^n(1 - u^n),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is a *linear* algebraic\n",
    "equation for the unknown value $u^{n+1}$ that we can easily solve:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+1} = u^n + \\Delta t\\,u^n(1 - u^n){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nonlinearity in the original equation poses in this case no difficulty\n",
    "in the discrete algebraic equation.\n",
    "Any other explicit scheme in time will also give only linear\n",
    "algebraic equations\n",
    "to solve. For example, a typical 2nd-order Runge-Kutta method\n",
    "for ([1](#nonlin:timediscrete:logistic:eq)) leads to the following\n",
    "formulas:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "u^* &= u^n + \\Delta t u^n(1 - u^n),\\\\\n",
    "u^{n+1} &= u^n + \\Delta t \\frac{1}{2} \\left(\n",
    "u^n(1 - u^n) + u^*(1 - u^*))\n",
    "\\right){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is linear in the unknown $u^*$. Then $u^*$ is\n",
    "known in the next step, which is linear in the unknown $u^{n+1}$ .\n",
    "For this scheme we have an explicit formula for the unknown \n",
    "and the scheme is therefore called an *explicit scheme*. \n",
    "\n",
    "\n",
    "## Exact solution of nonlinear algebraic equations\n",
    "<div id=\"nonlin:timediscrete:logistic:roots\"></div>\n",
    "\n",
    "Switching to a Backward Euler scheme for\n",
    "([1](#nonlin:timediscrete:logistic:eq)),"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:eq:BE\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n} - u^{n-1}}{\\Delta t} = u^n(1 - u^n),\n",
    "\\label{nonlin:timediscrete:logistic:eq:BE} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "results in a nonlinear algebraic equation for the unknown value $u^n$.\n",
    "The equation is of quadratic type:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Delta t (u^n)^2 + (1-\\Delta t)u^n - u^{n-1} = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and may be solved exactly by the well-known formula for such equations.\n",
    "Before we do so, however, we will\n",
    "introduce a shorter, and often cleaner, notation for\n",
    "nonlinear algebraic equations at a given time level. The notation is\n",
    "inspired by the natural notation (i.e., variable names) used in a\n",
    "program, especially in more advanced partial differential equation\n",
    "problems. The unknown in the algebraic equation is denoted by $u$,\n",
    "while $u^{(1)}$ is the value of the unknown at the previous time level\n",
    "(in general, $u^{(\\ell)}$ is the value of the unknown $\\ell$ levels\n",
    "back in time). The notation will be frequently used in later\n",
    "sections. What is meant by $u$ should be evident from the context: $u$\n",
    "may be 1) the exact solution of the ODE/PDE problem,\n",
    "2) the numerical approximation to the exact solution, or 3) the unknown\n",
    "solution at a certain time level.\n",
    "\n",
    "The quadratic equation for the unknown $u^n$ in\n",
    "([2](#nonlin:timediscrete:logistic:eq:BE)) can, with the new\n",
    "notation, be written"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:eq:F\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F(u) = \\Delta t u^2 + (1-\\Delta t)u - u^{(1)} = 0{\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:eq:F} \\tag{3}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is readily found to be"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:eq:roots\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u = \\frac{1}{2\\Delta t}\n",
    "\\left(-1+\\Delta t \\pm \\sqrt{(1-\\Delta t)^2 - 4\\Delta t u^{(1)}}\\right)\n",
    "{\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:eq:roots} \\tag{4}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we encounter a fundamental challenge with nonlinear\n",
    "algebraic equations:\n",
    "the equation may have more than one solution. How do we pick the right\n",
    "solution? This is in general a hard problem.\n",
    "In the present simple case, however, we can analyze the roots mathematically\n",
    "and provide an answer. The idea is to expand the roots\n",
    "in a series in $\\Delta t$ and truncate after the linear term since\n",
    "the Backward Euler scheme will introduce an error proportional to\n",
    "$\\Delta t$ anyway. Using `sympy` we find the following Taylor series\n",
    "expansions of the roots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "dt, u_1, u = sym.symbols('dt u_1 u')\n",
    "r1, r2 = sym.solve(dt*u**2 + (1-dt)*u - u_1, u)  # find roots\n",
    "r1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((r1.series(dt, 0, 2)))    # 2 terms in dt, around dt=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((r2.series(dt, 0, 2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the `r1` root, corresponding to\n",
    "a minus sign in front of the square root in\n",
    "([4](#nonlin:timediscrete:logistic:eq:roots)),\n",
    "behaves as $1/\\Delta t$ and will therefore\n",
    "blow up as $\\Delta t\\rightarrow 0$! Since we know that $u$ takes on\n",
    "finite values, actually it is less than or equal to 1,\n",
    "only the `r2` root is of relevance in this case: as $\\Delta t\\rightarrow 0$,\n",
    "$u\\rightarrow u^{(1)}$, which is the expected result.\n",
    "\n",
    "For those who are not well experienced with approximating mathematical\n",
    "formulas by series expansion, an alternative method of investigation\n",
    "is simply to compute the limits of the two roots as $\\Delta t\\rightarrow 0$\n",
    "and see if a limit unreasonable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((r1.limit(dt, 0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((r2.limit(dt, 0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linearization\n",
    "\n",
    "When the time integration of an ODE results in a nonlinear algebraic\n",
    "equation, we must normally find its solution by defining a sequence\n",
    "of linear equations and hope that the solutions of these linear equations\n",
    "converge to the desired solution of the nonlinear algebraic equation.\n",
    "Usually, this means solving the linear equation repeatedly in an\n",
    "iterative fashion.\n",
    "Alternatively, the nonlinear equation can sometimes be approximated by one\n",
    "linear equation, and consequently there is no need for iteration.\n",
    "\n",
    "\n",
    "Constructing a linear equation from a nonlinear one requires\n",
    "*linearization* of each nonlinear term. This can be done manually\n",
    "as in Picard iteration, or fully algorithmically as in Newton's method.\n",
    "Examples will best illustrate how to linearize nonlinear problems.\n",
    "\n",
    "\n",
    "## Picard iteration\n",
    "<div id=\"nonlin:timediscrete:logistic:Picard\"></div>\n",
    "\n",
    "\n",
    "Let us write ([3](#nonlin:timediscrete:logistic:eq:F)) in a\n",
    "more compact form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u) = au^2 + bu + c = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $a=\\Delta t$, $b=1-\\Delta t$, and $c=-u^{(1)}$.\n",
    "Let $u^{-}$ be an available approximation of the unknown $u$.\n",
    "\n",
    "Then we can linearize the term $u^2$ simply by writing\n",
    "$u^{-}u$. The resulting equation, $\\hat F(u)=0$, is now linear\n",
    "and hence easy to solve:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u)\\approx\\hat F(u) = au^{-}u + bu + c = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the equation $\\hat F=0$ is only approximate, the solution $u$\n",
    "does not equal the exact solution ${u_{\\small\\mbox{e}}}$ of the exact\n",
    "equation $F({u_{\\small\\mbox{e}}})=0$, but we can hope that $u$ is closer to\n",
    "${u_{\\small\\mbox{e}}}$ than $u^{-}$ is, and hence it makes sense to repeat the\n",
    "procedure, i.e., set $u^{-}=u$ and solve $\\hat F(u)=0$ again.\n",
    "There is no guarantee that $u$ is closer to ${u_{\\small\\mbox{e}}}$ than $u^{-}$,\n",
    "but this approach has proven to be effective in a wide range of\n",
    "applications.\n",
    "\n",
    "The idea of turning a nonlinear equation into a linear one by\n",
    "using an approximation $u^{-}$ of $u$ in the nonlinear terms is\n",
    "a widely used approach that goes under many names:\n",
    "*fixed-point iteration*, the method of *successive substitutions*,\n",
    "*nonlinear Richardson iteration*, and *Picard iteration*.\n",
    "We will stick to the latter name.\n",
    "\n",
    "\n",
    "Picard iteration for solving the nonlinear equation\n",
    "arising from the Backward Euler discretization of the logistic\n",
    "equation can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = -\\frac{c}{au^{-} + b},\\quad u^{-}\\ \\leftarrow\\ u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\leftarrow$ symbols means assignment (we set $u^{-}$ equal to\n",
    "the value of $u$).\n",
    "The iteration is started with the value of the unknown at the\n",
    "previous time level: $u^{-}=u^{(1)}$.\n",
    "\n",
    "\n",
    "Some prefer an explicit iteration counter as superscript\n",
    "in the mathematical notation. Let $u^k$ be the computed approximation\n",
    "to the solution in iteration $k$. In iteration $k+1$ we want\n",
    "to solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "au^k u^{k+1} + bu^{k+1} + c = 0\\quad\\Rightarrow\\quad u^{k+1}\n",
    "= -\\frac{c}{au^k + b},\\quad k=0,1,\\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we need to perform the iteration at every time level, the\n",
    "time level counter is often also included:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "au^{n,k} u^{n,k+1} + bu^{n,k+1} - u^{n-1} = 0\\quad\\Rightarrow\\quad u^{n,k+1}\n",
    "= \\frac{u^{n-1}}{au^{n,k} + b},\\quad k=0,1,\\ldots,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with the start value $u^{n,0}=u^{n-1}$ and the final converged value\n",
    "$u^{n}=u^{n,k}$ for sufficiently large $k$.\n",
    "\n",
    "However, we will normally apply a mathematical notation in our\n",
    "final formulas that is as close as possible to what we aim to write\n",
    "in a computer code and then it becomes natural to use $u$ and $u^{-}$\n",
    "instead of $u^{k+1}$ and $u^k$ or $u^{n,k+1}$ and $u^{n,k}$.\n",
    "\n",
    "\n",
    "\n",
    "### Stopping criteria\n",
    "\n",
    "The iteration method can typically be terminated when the change\n",
    "in the solution is smaller than a tolerance $\\epsilon_u$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "|u - u^{-}| \\leq\\epsilon_u,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or when the residual in the equation is sufficiently small ($< \\epsilon_r$),"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "|F(u)|= |au^2+bu + c| < \\epsilon_r{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A single Picard iteration\n",
    "\n",
    "Instead of iterating until a stopping criterion is fulfilled, one may\n",
    "iterate a specific number of times. Just one Picard iteration is\n",
    "popular as this corresponds to the intuitive idea of approximating a\n",
    "nonlinear term like $(u^n)^2$ by $u^{n-1}u^n$. This follows from the\n",
    "linearization $u^{-}u^n$ and the initial choice of $u^{-}=u^{n-1}$ at\n",
    "time level $t_n$. In other words, a single Picard iteration\n",
    "corresponds to using the solution at the previous time level to\n",
    "linearize nonlinear terms. The resulting discretization becomes (using\n",
    "proper values for $a$, $b$, and $c$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:BE:Picard:1it\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n} - u^{n-1}}{\\Delta t} = u^n(1 - u^{n-1}),\n",
    "\\label{nonlin:timediscrete:logistic:BE:Picard:1it} \\tag{5}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is a linear algebraic equation in the unknown $u^n$, making\n",
    "it easy to solve for $u^n$ without any need for\n",
    "any alternative notation.\n",
    "\n",
    "We shall later refer to the strategy of taking one Picard step, or\n",
    "equivalently, linearizing terms with use of the solution at the\n",
    "previous time step, as the *Picard1* method. It is a widely used\n",
    "approach in science and technology, but with some limitations if\n",
    "$\\Delta t$ is not sufficiently small (as will be illustrated later).\n",
    "\n",
    "\n",
    "**Notice.**\n",
    "\n",
    "\n",
    "Equation ([5](#nonlin:timediscrete:logistic:BE:Picard:1it)) does not\n",
    "correspond to a \"pure\" finite difference method where the equation\n",
    "is sampled at a point and derivatives replaced by differences (because\n",
    "the $u^{n-1}$ term on the right-hand side must then be $u^n$). The\n",
    "best interpretation of the scheme\n",
    "([5](#nonlin:timediscrete:logistic:BE:Picard:1it)) is a Backward Euler\n",
    "difference combined with a single (perhaps insufficient) Picard\n",
    "iteration at each time level, with the value at the previous time\n",
    "level as start for the Picard iteration.\n",
    "\n",
    "\n",
    "\n",
    "## Linearization by a geometric mean\n",
    "<div id=\"nonlin:timediscrete:logistic:geometric:mean\"></div>\n",
    "\n",
    "We consider now a Crank-Nicolson discretization of\n",
    "([1](#nonlin:timediscrete:logistic:eq)). This means that the\n",
    "time derivative is approximated by a centered\n",
    "difference,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_t u = u(1-u)]^{n+\\frac{1}{2}},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "written out as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:geometric:mean:scheme\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n+1}-u^n}{\\Delta t} = u^{n+\\frac{1}{2}} -\n",
    "(u^{n+\\frac{1}{2}})^2{\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:geometric:mean:scheme} \\tag{6}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first term $u^{n+\\frac{1}{2}}$ is normally approximated by an arithmetic\n",
    "mean,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+\\frac{1}{2}}\\approx \\frac{1}{2}(u^n + u^{n+1}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "such that the scheme involves the unknown function only at the time levels\n",
    "where we actually compute it.\n",
    "The same arithmetic mean applied to the second term gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u^{n+\\frac{1}{2}})^2\\approx \\frac{1}{4}(u^n + u^{n+1})^2,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is nonlinear in the unknown $u^{n+1}$.\n",
    "However, using a *geometric mean* for $(u^{n+\\frac{1}{2}})^2$\n",
    "is a way of linearizing the nonlinear term in\n",
    "([6](#nonlin:timediscrete:logistic:geometric:mean:scheme)):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u^{n+\\frac{1}{2}})^2\\approx u^nu^{n+1}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using an arithmetic mean on the linear $u^{n+\\frac{1}{2}}$ term in\n",
    "([6](#nonlin:timediscrete:logistic:geometric:mean:scheme)) and a geometric\n",
    "mean for the second term, results in a linearized equation for the\n",
    "unknown $u^{n+1}$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n+1}-u^n}{\\Delta t} =\n",
    "\\frac{1}{2}(u^n + u^{n+1}) - u^nu^{n+1},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can readily be solved:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+1} = \\frac{1 + \\frac{1}{2}\\Delta t}{1+\\Delta t u^n - \\frac{1}{2}\\Delta t}\n",
    "u^n{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This scheme can be coded directly, and since\n",
    "there is no nonlinear algebraic equation to iterate over,\n",
    "we skip the simplified notation with $u$ for $u^{n+1}$\n",
    "and $u^{(1)}$ for $u^n$. The technique with using\n",
    "a geometric average is an example of transforming a nonlinear\n",
    "algebraic equation to a linear one, without any need for iterations.\n",
    "\n",
    "The geometric mean approximation is often very effective for\n",
    "linearizing quadratic nonlinearities. Both the arithmetic and geometric mean\n",
    "approximations have truncation errors of order $\\Delta t^2$ and are\n",
    "therefore compatible with the truncation error ${\\mathcal{O}(\\Delta t^2)}$\n",
    "of the centered difference approximation for $u^\\prime$ in the Crank-Nicolson\n",
    "method.\n",
    "\n",
    "Applying the operator notation for the means and finite differences,\n",
    "the linearized Crank-Nicolson scheme for the logistic equation can be\n",
    "compactly expressed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_t u = \\overline{u}^{t} + \\overline{u^2}^{t,g}]^{n+\\frac{1}{2}}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remark.**\n",
    "\n",
    "If we use an arithmetic instead of a geometric mean\n",
    "for the nonlinear term in\n",
    "([6](#nonlin:timediscrete:logistic:geometric:mean:scheme)),\n",
    "we end up with a nonlinear term $(u^{n+1})^2$.\n",
    "This term can be linearized as $u^{-}u^{n+1}$ in a Picard iteration\n",
    "approach and in particular as\n",
    "$u^nu^{n+1}$ in a Picard1 iteration approach.\n",
    "The latter gives a scheme almost identical to the one arising from\n",
    "a geometric mean (the difference in $u^{n+1}$\n",
    "being $\\frac{1}{4}\\Delta t u^n(u^{n+1}-u^n)\\approx \\frac{1}{4}\\Delta t^2\n",
    "u^\\prime u$, i.e., a difference of size $\\Delta t^2$).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Newton's method\n",
    "<div id=\"nonlin:timediscrete:logistic:Newton\"></div>\n",
    "\n",
    "\n",
    "The Backward Euler scheme ([2](#nonlin:timediscrete:logistic:eq:BE))\n",
    "for the logistic equation leads to a nonlinear algebraic equation\n",
    "([3](#nonlin:timediscrete:logistic:eq:F)). Now we write any nonlinear\n",
    "algebraic equation in the general and compact form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u) = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton's method linearizes this equation by approximating $F(u)$ by\n",
    "its Taylor series expansion around a computed value $u^{-}$\n",
    "and keeping only the linear part:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "F(u) &= F(u^{-}) + F^{\\prime}(u^{-})(u - u^{-}) + {\\frac{1}{2}}F^{\\prime\\prime}(u^{-})(u-u^{-})^2\n",
    "+\\cdots\\\\\n",
    "& \\approx F(u^{-}) + F^{\\prime}(u^{-})(u - u^{-}) = \\hat F(u){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear equation $\\hat F(u)=0$ has the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = u^{-} - \\frac{F(u^{-})}{F^{\\prime}(u^{-})}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expressed with an iteration index in the unknown, Newton's method takes\n",
    "on the more familiar mathematical form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{k+1} = u^k - \\frac{F(u^k)}{F^{\\prime}(u^k)},\\quad k=0,1,\\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the method converges, it can be shown that the error in iteration $k+1$ of Newton's method is\n",
    "proportional to\n",
    "the square of the error in iteration $k$, a result referred to as\n",
    "*quadratic convergence*. This means that for\n",
    "small errors the method converges very fast, and in particular much\n",
    "faster than Picard iteration and other iteration methods.\n",
    "(The proof of this result is found in most textbooks on numerical analysis.)\n",
    "However, the quadratic convergence appears only if $u^k$ is sufficiently\n",
    "close to the solution. Further away from the solution the method can\n",
    "easily converge very slowly or diverge. The reader is encouraged to do\n",
    "[nonlin:exer:Newton:problems1](#nonlin:exer:Newton:problems1) to get a better understanding\n",
    "for the behavior of the method.\n",
    "\n",
    "Application of Newton's method to the logistic equation discretized\n",
    "by the Backward Euler method is straightforward\n",
    "as we have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u) = au^2 + bu + c,\\quad a=\\Delta t,\\ b = 1-\\Delta t,\\ c=-u^{(1)},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F^{\\prime}(u) = 2au + b{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The iteration method becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:Newton:alg1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u = u^{-} + \\frac{a(u^{-})^2 + bu^{-} + c}{2au^{-} + b},\\quad\n",
    "u^{-}\\ \\leftarrow u{\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:Newton:alg1} \\tag{7}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At each time level, we start the iteration by setting $u^{-}=u^{(1)}$.\n",
    "Stopping criteria as listed for the Picard iteration can be used also\n",
    "for Newton's method.\n",
    "\n",
    "An alternative mathematical form, where we write out $a$, $b$, and $c$,\n",
    "and use a time level counter $n$ and an iteration counter $k$, takes\n",
    "the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:Newton:alg2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u^{n,k+1} = u^{n,k} +\n",
    "\\frac{\\Delta t (u^{n,k})^2 + (1-\\Delta t)u^{n,k} - u^{n-1}}\n",
    "{2\\Delta t u^{n,k} + 1 - \\Delta t},\\quad u^{n,0}=u^{n-1},\n",
    "\\label{nonlin:timediscrete:logistic:Newton:alg2} \\tag{8}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $k=0,1,\\ldots$.\n",
    "A program implementation is much closer to ([7](#nonlin:timediscrete:logistic:Newton:alg1)) than to ([8](#nonlin:timediscrete:logistic:Newton:alg2)), but\n",
    "the latter is better aligned with the established mathematical\n",
    "notation used in the literature.\n",
    "\n",
    "## Relaxation\n",
    "<div id=\"nonlin:timediscrete:logistic:relaxation\"></div>\n",
    "\n",
    "\n",
    "One iteration in Newton's method or\n",
    "Picard iteration consists of solving a linear problem $\\hat F(u)=0$.\n",
    "Sometimes convergence problems arise because the new solution $u$\n",
    "of $\\hat F(u)=0$ is \"too far away\" from the previously computed\n",
    "solution $u^{-}$. A remedy is to introduce a relaxation, meaning that\n",
    "we first solve $\\hat F(u^*)=0$ for a suggested value $u^*$ and\n",
    "then we take $u$ as a weighted mean of what we had, $u^{-}$, and\n",
    "what our linearized equation $\\hat F=0$ suggests, $u^*$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = \\omega u^* + (1-\\omega) u^{-}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter $\\omega$\n",
    "is known as a *relaxation parameter*, and a choice $\\omega < 1$\n",
    "may prevent divergent iterations.\n",
    "\n",
    "Relaxation in Newton's method can be directly incorporated\n",
    "in the basic iteration formula:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:timediscrete:logistic:relaxation:Newton:formula\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u = u^{-} - \\omega \\frac{F(u^{-})}{F^{\\prime}(u^{-})}{\\thinspace .}\n",
    "\\label{nonlin:timediscrete:logistic:relaxation:Newton:formula} \\tag{9}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation and experiments\n",
    "<div id=\"nonlin:timediscrete:logistic:impl\"></div>\n",
    "\n",
    "\n",
    "The program `logistic.py` contains\n",
    "implementations of all the methods described above.\n",
    "Below is an extract of the file showing how the Picard and Newton\n",
    "methods are implemented for a Backward Euler discretization of\n",
    "the logistic equation.\n",
    "\n",
    "<!-- @@@CODE src/logistic.py fromto: def BE_logistic@def CN_logistic -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BE_logistic(u0, dt, Nt, choice='Picard',\n",
    "                eps_r=1E-3, omega=1, max_iter=1000):\n",
    "    if choice == 'Picard1':\n",
    "        choice = 'Picard'\n",
    "        max_iter = 1\n",
    "\n",
    "    u = np.zeros(Nt+1)\n",
    "    iterations = []\n",
    "    u[0] = u0\n",
    "    for n in range(1, Nt+1):\n",
    "        a = dt\n",
    "        b = 1 - dt\n",
    "        c = -u[n-1]\n",
    "\n",
    "        if choice == 'Picard':\n",
    "\n",
    "            def F(u):\n",
    "                return a*u**2 + b*u + c\n",
    "\n",
    "            u_ = u[n-1]\n",
    "            k = 0\n",
    "            while abs(F(u_)) > eps_r and k < max_iter:\n",
    "                u_ = omega*(-c/(a*u_ + b)) + (1-omega)*u_\n",
    "                k += 1\n",
    "            u[n] = u_\n",
    "            iterations.append(k)\n",
    "\n",
    "        elif choice == 'Newton':\n",
    "\n",
    "            def F(u):\n",
    "                return a*u**2 + b*u + c\n",
    "\n",
    "            def dF(u):\n",
    "                return 2*a*u + b\n",
    "\n",
    "            u_ = u[n-1]\n",
    "            k = 0\n",
    "            while abs(F(u_)) > eps_r and k < max_iter:\n",
    "                u_ = u_ - F(u_)/dF(u_)\n",
    "                k += 1\n",
    "            u[n] = u_\n",
    "            iterations.append(k)\n",
    "    return u, iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Crank-Nicolson method utilizing a linearization based on the\n",
    "geometric mean gives a simpler algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CN_logistic(u0, dt, Nt):\n",
    "    u = np.zeros(Nt+1)\n",
    "    u[0] = u0\n",
    "    for n in range(0, Nt):\n",
    "        u[n+1] = (1 + 0.5*dt)/(1 + dt*u[n] - 0.5*dt)*u[n]\n",
    "    return u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may run experiments with the model problem\n",
    "([1](#nonlin:timediscrete:logistic:eq)) and the different strategies for\n",
    "dealing with nonlinearities as described above. For a quite coarse\n",
    "time resolution, $\\Delta t=0.9$, use of a tolerance $\\epsilon_r=0.05$\n",
    "in the stopping criterion introduces an iteration error, especially in\n",
    "the Picard iterations, that is visibly much larger than the\n",
    "time discretization error due to a large $\\Delta t$. This is illustrated\n",
    "by comparing the upper two plots in\n",
    "[Figure](#nonlin:timediscrete:logistic:impl:fig:u). The one to\n",
    "the right has a stricter tolerance $\\epsilon = 10^{-3}$, which leads\n",
    "to all the curves involving backward Euler, using iterative solution by\n",
    "either Picard or Newton iterations, to be\n",
    "on top of each other (and no changes can be visually observed by\n",
    "reducing $\\epsilon_r$ further). The reason why Newton's method does\n",
    "much better than Picard iteration in the upper left plot is that\n",
    "Newton's method with one step comes far below the $\\epsilon_r$ tolerance,\n",
    "while the Picard iteration needs on average 7 iterations to bring the\n",
    "residual down to $\\epsilon_r=10^{-1}$, which gives insufficient\n",
    "accuracy in the solution of the nonlinear equation. It is obvious\n",
    "that the Picard1 method gives significant errors in addition to\n",
    "the time discretization unless the time step is as small as in\n",
    "the lower right plot.\n",
    "\n",
    "The *BE exact* curve corresponds to using the exact solution of the\n",
    "quadratic equation at each time level, so this curve is only affected\n",
    "by the Backward Euler time discretization.  The *CN gm* curve\n",
    "corresponds to the theoretically more accurate Crank-Nicolson\n",
    "discretization, combined with a geometric mean for linearization.\n",
    "Visually, this curve appears more accurate in all the plots, especially if we take the plot in\n",
    "the lower right with a small $\\Delta t$ and an appropriately small\n",
    "$\\epsilon_r$ value as the reference curve.\n",
    "\n",
    "\n",
    "When it comes to the need for iterations, [Figure](#nonlin:timediscrete:logistic:impl:fig:iter) displays the number of\n",
    "iterations required at each time level for Newton's method and\n",
    "Picard iteration. The smaller $\\Delta t$ is, the better starting value\n",
    "we have for the iteration, and the faster the convergence is.\n",
    "With $\\Delta t = 0.9$ Picard iteration requires on average 32 iterations\n",
    "per time step for the stricter convergence criterion, but this number is dramatically reduced as $\\Delta t$\n",
    "is reduced.\n",
    "\n",
    "However, introducing relaxation and a parameter $\\omega=0.8$\n",
    "immediately reduces the average of 32 to 7, indicating that for the large\n",
    "$\\Delta t=0.9$, Picard iteration takes too long steps. An approximately optimal\n",
    "value for $\\omega$ in this case is 0.5, which results in an average of only\n",
    "2 iterations! An even more dramatic impact of $\\omega$ appears when\n",
    "$\\Delta t = 1$: Picard iteration does not convergence in 1000 iterations,\n",
    "but $\\omega=0.5$ again brings the average number of iterations down to 2.\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/logistic_u.png, width=800 frac=1] Impact of solution strategy and time step length on the solution. <div id=\"nonlin:timediscrete:logistic:impl:fig:u\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"nonlin:timediscrete:logistic:impl:fig:u\"></div>\n",
    "\n",
    "<p>Impact of solution strategy and time step length on the solution.</p>\n",
    "<img src=\"fig/logistic_u.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/logistic_iter.png, width=800 frac=1] Comparison of the number of iterations at various time levels for Picard and Newton iteration. <div id=\"nonlin:timediscrete:logistic:impl:fig:iter\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"nonlin:timediscrete:logistic:impl:fig:iter\"></div>\n",
    "\n",
    "<p>Comparison of the number of iterations at various time levels for Picard and Newton iteration.</p>\n",
    "<img src=\"fig/logistic_iter.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Generalization to a general nonlinear ODE\n",
    "<div id=\"nonlin:ode:generic\"></div>\n",
    "\n",
    "Let us see how the various methods in the previous sections\n",
    "can be applied to the more generic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:ode:generic:model\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u^{\\prime} = f(u, t),\n",
    "\\label{nonlin:ode:generic:model} \\tag{10}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $f$ is a nonlinear function of $u$.\n",
    "\n",
    "### Explicit time discretization\n",
    "\n",
    "Explicit ODE methods like the Forward Euler scheme, various Runge-Kutta methods,\n",
    "Adams-Bashforth methods all evaluate $f$ at time levels where\n",
    "$u$ is already computed, so nonlinearities in $f$ do not\n",
    "pose any difficulties.\n",
    "\n",
    "### Backward Euler discretization\n",
    "\n",
    "Approximating $u^{\\prime}$ by a backward difference leads to a Backward Euler\n",
    "scheme, which can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u^n) = u^{n} - \\Delta t\\, f(u^n, t_n) - u^{n-1}=0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or alternatively"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u) = u - \\Delta t\\, f(u, t_n) - u^{(1)} = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple Picard iteration, not knowing anything about the nonlinear\n",
    "structure of $f$, must approximate $f(u,t_n)$ by $f(u^{-},t_n)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\hat F(u) = u - \\Delta t\\, f(u^{-},t_n) - u^{(1)}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The iteration starts with $u^{-}=u^{(1)}$ and proceeds with repeating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^* = \\Delta t\\, f(u^{-},t_n) + u^{(1)},\\quad u = \\omega u^* + (1-\\omega)u^{-},\n",
    "\\quad u^{-}\\ \\leftarrow\\ u,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "until a stopping criterion is fulfilled.\n",
    "\n",
    "**Explicit vs implicit treatment of nonlinear terms.**\n",
    "\n",
    "Evaluating $f$ for a known $u^{-}$ is referred to as *explicit* treatment of\n",
    "$f$, while if $f(u,t)$ has some structure, say $f(u,t) = u^3$, parts of\n",
    "$f$ can involve the unknown $u$, as in the manual linearization\n",
    "like $(u^{-})^2u$, and then the treatment of $f$ is \"more implicit\"\n",
    "and \"less explicit\". This terminology is inspired by time discretization\n",
    "of $u^{\\prime}=f(u,t)$, where evaluating $f$ for known $u$ values gives\n",
    "explicit formulas for the unknown and hence explicit schemes, while treating $f$ \n",
    "or parts of $f$ implicitly, meaning that equations must be solved in terms of the \n",
    "unknown, makes $f$ contribute to the unknown terms in the equation at the new\n",
    "time level.\n",
    "\n",
    "Explicit treatment of $f$ usually means stricter conditions on\n",
    "$\\Delta t$ to achieve stability of time discretization schemes. The same\n",
    "applies to iteration techniques for nonlinear algebraic equations: the \"less\"\n",
    "we linearize $f$ (i.e., the more we keep of $u$ in the original formula),\n",
    "the faster the convergence may be.\n",
    "\n",
    "We may say that $f(u,t)=u^3$ is treated explicitly if we evaluate $f$\n",
    "as $(u^{-})^3$, semi or partially implicit if we linearize as $(u^{-})^2u$\n",
    "and fully implicit if we represent $f$ by $u^3$. (Of course, the\n",
    "fully implicit representation will require further linearization,\n",
    "but with $f(u,t)=u^2$ a fully implicit treatment is possible if\n",
    "the resulting quadratic equation is solved with a formula.)\n",
    "\n",
    "For the ODE $u^{\\prime}=-u^3$ with $f(u,t)=-u^3$ and coarse\n",
    "time resolution $\\Delta t = 0.4$, Picard iteration with $(u^{-})^2u$\n",
    "requires 8 iterations with $\\epsilon_r = 10^{-3}$ for the first\n",
    "time step, while $(u^{-})^3$ leads to 22 iterations. After about 10\n",
    "time steps both approaches are down to about 2 iterations per time\n",
    "step, but this example shows a potential of treating $f$ more\n",
    "implicitly.\n",
    "\n",
    "A trick to treat $f$ implicitly in Picard iteration is to\n",
    "evaluate it as $f(u^{-},t)u/u^{-}$. For a polynomial $f$, $f(u,t)=u^m$,\n",
    "this corresponds to $(u^{-})^{m}u/u^{-}=(u^{-})^{m-1}u$. Sometimes this more implicit\n",
    "treatment has no effect, as with $f(u,t)=\\exp(-u)$ and $f(u,t)=\\ln (1+u)$,\n",
    "but with $f(u,t)=\\sin(2(u+1))$, the $f(u^{-},t)u/u^{-}$ trick\n",
    "leads to 7, 9, and 11 iterations during the first three steps, while\n",
    "$f(u^{-},t)$ demands 17, 21, and 20 iterations.\n",
    "(Experiments can be done with the code [`ODE_Picard_tricks.py`](${fem_src}/ODE_Picard_tricks.py).)\n",
    "\n",
    "\n",
    "\n",
    "Newton's method applied to a Backward Euler discretization of\n",
    "$u^{\\prime}=f(u,t)$\n",
    "requires the computation of the derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F^{\\prime}(u) = 1 - \\Delta t\\frac{\\partial f}{\\partial u}(u,t_n){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting with the solution at the previous time level, $u^{-}=u^{(1)}$,\n",
    "we can just use the standard formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:ode:generic:Newton\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u = u^{-} - \\omega \\frac{F(u^{-})}{F^{\\prime}(u^{-})}\n",
    "= u^{-} - \\omega \\frac{u^{-} - \\Delta t\\, f(u^{-}, t_n) - u^{(1)}}{1 - \\Delta t\n",
    "\\frac{\\partial}{\\partial u}f(u^{-},t_n)}\n",
    "{\\thinspace .}\n",
    "\\label{nonlin:ode:generic:Newton} \\tag{11}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- The geometric mean trick cannot be used unless we know that $f$ has -->\n",
    "<!-- a special structure with quadratic expressions in $u$. -->\n",
    "\n",
    "### Crank-Nicolson discretization\n",
    "\n",
    "The standard Crank-Nicolson scheme with arithmetic mean approximation of\n",
    "$f$ takes the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n+1} - u^n}{\\Delta t} = \\frac{1}{2}(f(u^{n+1}, t_{n+1})\n",
    "+ f(u^n, t_n)){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write the scheme as a nonlinear algebraic equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:ode:generic:Newton2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F(u) = u - u^{(1)} - \\Delta t{\\frac{1}{2}}f(u,t_{n+1}) -\n",
    "\\Delta t{\\frac{1}{2}}f(u^{(1)},t_{n}) = 0{\\thinspace .}\n",
    "\\label{nonlin:ode:generic:Newton2} \\tag{12}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Picard iteration scheme must in general employ the linearization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\hat F(u) = u - u^{(1)} - \\Delta t{\\frac{1}{2}}f(u^{-},t_{n+1}) -\n",
    "\\Delta t{\\frac{1}{2}}f(u^{(1)},t_{n}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "while Newton's method can apply the general formula\n",
    "([11](#nonlin:ode:generic:Newton)) with $F(u)$ given in\n",
    "([12](#nonlin:ode:generic:Newton2)) and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F^{\\prime}(u)= 1 - \\frac{1}{2}\\Delta t\\frac{\\partial f}{\\partial u}(u,t_{n+1}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- What about pendulum sin(u) as u/u_ sin(u_)? Check in odespy if it -->\n",
    "<!-- converges faster (should be able to store the no of Newton and -->\n",
    "<!-- Picard iterations in the classes and poll afterwards). It the trick -->\n",
    "<!-- pays off, describe it here. Can odespy be used here? That is, can we -->\n",
    "<!-- provide the linearization? No...? -->\n",
    "\n",
    "## Systems of ODEs\n",
    "<div id=\"nonlin:ode:generic:sys:pendulum\"></div>\n",
    "\n",
    "We may write a system of ODEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{d}{dt}u_0(t) &= f_0(u_0(t),u_1(t),\\ldots,u_N(t),t),\\\\\n",
    "\\frac{d}{dt}u_1(t) &= f_1(u_0(t),u_1(t),\\ldots,u_N(t),t),\\\\\n",
    "&\\vdots\\\\\n",
    "\\frac{d}{dt}u_m(t) &= f_m(u_0(t),u_1(t),\\ldots,u_N(t),t),\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u^{\\prime} = f(u,t),\\quad u(0)=U_0,\n",
    "\\label{_auto1} \\tag{13}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "if we interpret $u$ as a vector $u=(u_0(t),u_1(t),\\ldots,u_N(t))$\n",
    "and $f$ as a vector function with components\n",
    "$(f_0(u,t),f_1(u,t),\\ldots,f_N(u,t))$.\n",
    "\n",
    "\n",
    "\n",
    "Most solution methods for scalar ODEs, including\n",
    "the Forward and Backward Euler schemes and the\n",
    "Crank-Nicolson method, generalize in a\n",
    "straightforward way to systems of ODEs simply by using vector\n",
    "arithmetics instead of scalar arithmetics, which corresponds to\n",
    "applying the scalar scheme to each component of the system.  For\n",
    "example, here is a backward difference scheme applied to each\n",
    "component,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{u_0^n- u_0^{n-1}}{\\Delta t} &= f_0(u^n,t_n),\\\\\n",
    "\\frac{u_1^n- u_1^{n-1}}{\\Delta t} &= f_1(u^n,t_n),\\\\\n",
    "&\\vdots\\\\\n",
    "\\frac{u_N^n- u_N^{n-1}}{\\Delta t} &= f_N(u^n,t_n),\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can be written more compactly in vector form as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^n- u^{n-1}}{\\Delta t} = f(u^n,t_n){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a *system of algebraic equations*,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^n - \\Delta t\\,f(u^n,t_n) - u^{n-1}=0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or written out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "u_0^n - \\Delta t\\, f_0(u^n,t_n) - u_0^{n-1} &= 0,\\\\\n",
    "&\\vdots\\\\\n",
    "u_N^n - \\Delta t\\, f_N(u^n,t_n) - u_N^{n-1} &= 0{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "We shall address the $2\\times 2$ ODE system for\n",
    "oscillations of a pendulum\n",
    "subject to gravity and air drag. The system can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\dot\\omega = -\\sin\\theta -\\beta \\omega |\\omega|,\n",
    "\\label{_auto2} \\tag{14}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\dot\\theta = \\omega,\n",
    "\\label{_auto3} \\tag{15}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\beta$ is a dimensionless parameter (this is the scaled, dimensionless\n",
    "version of the original, physical model). The unknown components of the\n",
    "system are the\n",
    "angle $\\theta(t)$ and the angular velocity $\\omega(t)$.\n",
    "We introduce $u_0=\\omega$ and $u_1=\\theta$, which leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "u_0^{\\prime} = f_0(u,t) &= -\\sin u_1 - \\beta u_0|u_0|,\\\\\n",
    "u_1^{\\prime} = f_1(u,t) &= u_0{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Crank-Nicolson scheme reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u_0^{n+1}-u_0^{n}}{\\Delta t} = -\\sin u_1^{n+\\frac{1}{2}}\n",
    "- \\beta u_0^{n+\\frac{1}{2}}|u_0^{n+\\frac{1}{2}}|\\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    " \\approx -\\sin\\left(\\frac{1}{2}(u_1^{n+1} + u_1^n)\\right)\n",
    "- \\beta\\frac{1}{4} (u_0^{n+1} + u_0^n)|u_0^{n+1}+u_0^n|,\n",
    "\\label{_auto4} \\tag{16}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto5\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\frac{u_1^{n+1}-u_1^n}{\\Delta t} = u_0^{n+\\frac{1}{2}}\\approx\n",
    "\\frac{1}{2} (u_0^{n+1}+u_0^n){\\thinspace .}\n",
    "\\label{_auto5} \\tag{17}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a *coupled system* of two nonlinear algebraic equations\n",
    "in two unknowns $u_0^{n+1}$ and $u_1^{n+1}$.\n",
    "\n",
    "Using the notation $u_0$ and $u_1$ for the unknowns $u_0^{n+1}$ and\n",
    "$u_1^{n+1}$ in this system, writing $u_0^{(1)}$ and\n",
    "$u_1^{(1)}$ for the previous values $u_0^n$ and $u_1^n$, multiplying\n",
    "by $\\Delta t$ and moving the terms to the left-hand sides, gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:ode:generic:sys:pendulum:u0\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u_0 - u_0^{(1)} + \\Delta t\\,\\sin\\left(\\frac{1}{2}(u_1 + u_1^{(1)})\\right)\n",
    "+ \\frac{1}{4}\\Delta t\\beta (u_0 + u_0^{(1)})|u_0 + u_0^{(1)}| =0,\n",
    "\\label{nonlin:ode:generic:sys:pendulum:u0} \\tag{18}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:ode:generic:sys:pendulum:u1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "u_1 - u_1^{(1)} -\\frac{1}{2}\\Delta t(u_0 + u_0^{(1)}) =0{\\thinspace .}\n",
    "\\label{nonlin:ode:generic:sys:pendulum:u1} \\tag{19}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously, we have a need for solving systems of nonlinear algebraic\n",
    "equations, which is the topic of the next section.\n",
    "\n",
    "# Systems of nonlinear algebraic equations\n",
    "<div id=\"nonlin:systems:alg\"></div>\n",
    "\n",
    "Implicit time discretization methods for a system of ODEs, or a PDE,\n",
    "lead to *systems* of nonlinear algebraic equations, written\n",
    "compactly as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u) = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $u$ is a vector of unknowns $u=(u_0,\\ldots,u_N)$, and\n",
    "$F$ is a vector function: $F=(F_0,\\ldots,F_N)$.\n",
    "The system at the end of the section [Systems of ODEs](#nonlin:ode:generic:sys:pendulum) fits\n",
    "this notation with $N=2$, $F_0(u)$ given by the left-hand side of\n",
    "([18](#nonlin:ode:generic:sys:pendulum:u0)), while $F_1(u)$ is\n",
    "the left-hand side of ([19](#nonlin:ode:generic:sys:pendulum:u1)).\n",
    "\n",
    "Sometimes the equation system has a special structure because of the\n",
    "underlying problem, e.g.,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A(u)u = b(u),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $A(u)$ as\n",
    "an $(N+1)\\times (N+1)$ matrix function of $u$ and $b$ as a vector\n",
    "function: $b=(b_0,\\ldots,b_N)$.\n",
    "\n",
    "We shall next explain how Picard iteration and Newton's method\n",
    "can be applied to systems like $F(u)=0$ and $A(u)u=b(u)$.\n",
    "The exposition has a focus on ideas and practical computations.\n",
    "More theoretical considerations, including quite general results\n",
    "on convergence properties\n",
    "of these methods, can be found in Kelley [[Kelley_1995]](#Kelley_1995).\n",
    "\n",
    "## Picard iteration\n",
    "<div id=\"nonlin:systems:alg:Picard\"></div>\n",
    "\n",
    "We cannot apply Picard iteration to nonlinear equations unless there is\n",
    "some special structure. For the commonly arising case\n",
    "$A(u)u=b(u)$ we can linearize the\n",
    "product $A(u)u$ to $A(u^{-})u$ and $b(u)$ as $b(u^{-})$.\n",
    "That is, we use the most previously\n",
    "computed approximation in $A$ and $b$ to arrive at a *linear system* for\n",
    "$u$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A(u^{-})u = b(u^{-}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A relaxed iteration takes the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A(u^{-})u^* = b(u^{-}),\\quad u = \\omega u^* + (1-\\omega)u^{-}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, we solve a system of nonlinear algebraic equations as\n",
    "a sequence of linear systems.\n",
    "\n",
    "**Algorithm for relaxed Picard iteration.**\n",
    "\n",
    "Given $A(u)u=b(u)$ and an initial guess $u^{-}$, iterate until convergence:\n",
    "\n",
    "1. solve $A(u^{-})u^* = b(u^{-})$ with respect to $u^*$\n",
    "\n",
    "2. $u = \\omega u^* + (1-\\omega) u^{-}$\n",
    "\n",
    "3. $u^{-}\\ \\leftarrow\\ u$\n",
    "\n",
    "\n",
    "\n",
    "\"Until convergence\" means that the iteration is stopped when the\n",
    "change in the unknown, $||u - u^{-}||$, or the residual $||A(u)u-b||$,\n",
    "is sufficiently small, see the section [Stopping criteria](#nonlin:systems:alg:terminate) for\n",
    "more details.\n",
    "\n",
    "## Newton's method\n",
    "<div id=\"nonlin:systems:alg:Newton\"></div>\n",
    "\n",
    "The natural starting point for Newton's method is the general\n",
    "nonlinear vector equation $F(u)=0$.\n",
    "As for a scalar equation, the idea is to approximate $F$\n",
    "around a known value $u^{-}$ by a linear function $\\hat F$,\n",
    "calculated from the first two terms of a Taylor expansion of\n",
    "$F$. In the multi-variate case these two terms become"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F(u^{-}) + J(u^{-}) \\cdot (u - u^{-}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $J$ is the *Jacobian* of $F$, defined by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "J_{i,j} = \\frac{\\partial F_i}{\\partial u_j}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the original nonlinear system is approximated by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\hat F(u) = F(u^{-}) + J(u^{-}) \\cdot (u - u^{-})=0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is linear in $u$ and can be solved in a two-step procedure:\n",
    "first solve $J\\delta u = -F(u^{-})$ with respect to the vector $\\delta u$\n",
    "and then update $u = u^{-} + \\delta u$.\n",
    "A relaxation parameter can easily be incorporated:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = \\omega(u^{-} +\\delta u)\n",
    "+ (1-\\omega)u^{-} = u^{-}  + \\omega\\delta u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Algorithm for Newton's method.**\n",
    "\n",
    "Given $F(u)=0$ and an initial guess $u^{-}$, iterate until convergence:\n",
    "\n",
    "1. solve $J\\delta u = -F(u^{-})$ with respect to $\\delta u$\n",
    "\n",
    "2. $u = u^{-} + \\omega\\delta u$\n",
    "\n",
    "3. $u^{-}\\ \\leftarrow\\ u$\n",
    "\n",
    "\n",
    "\n",
    "For the special system with structure $A(u)u=b(u)$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_i = \\sum_k A_{i,k}(u)u_k - b_i(u),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "one gets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto6\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "J_{i,j} = \\sum_k \\frac{\\partial A_{i,k}}{\\partial u_j}u_k\n",
    "+ A_{i,j} -\n",
    "\\frac{\\partial b_i}{\\partial u_j}{\\thinspace .}\n",
    "\\label{_auto6} \\tag{20}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We realize that the Jacobian needed in Newton's method consists of\n",
    "$A(u^{-})$ as in the Picard iteration plus two additional terms\n",
    "arising from the differentiation. Using the notation $A^{\\prime}(u)$ for\n",
    "$\\partial A/\\partial u$ (a quantity with three indices: $\\partial\n",
    "A_{i,k}/\\partial u_j$), and $b^{\\prime}(u)$ for $\\partial b/\\partial u$ (a\n",
    "quantity with two indices: $\\partial b_i/\\partial u_j$), we can write\n",
    "the linear system to be solved as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(A + A^{\\prime}u + b^{\\prime})\\delta u = -Au + b,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(A(u^{-}) + A^{\\prime}(u^{-})u^{-} + b^{\\prime}(u^{-}))\\delta u\n",
    "= -A(u^{-})u^{-} + b(u^{-}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rearranging the terms demonstrates the difference from the system\n",
    "solved in each Picard iteration:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{A(u^{-})(u^{-}+\\delta u) - b(u^{-})}_{\\hbox{Picard system}}\n",
    "+\\, \\gamma (A^{\\prime}(u^{-})u^{-} + b^{\\prime}(u^{-}))\\delta u\n",
    "= 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have inserted a parameter $\\gamma$ such that $\\gamma=0$\n",
    "gives the Picard system and $\\gamma=1$ gives the Newton system.\n",
    "Such a parameter can be handy in software to easily switch between\n",
    "the methods.\n",
    "\n",
    "**Combined algorithm for Picard and Newton iteration.**\n",
    "\n",
    "Given $A(u)$, $b(u)$, and an initial guess $u^{-}$, iterate until convergence:\n",
    "\n",
    "1. solve $(A + \\gamma(A^{\\prime}(u^{-})u^{-} +\n",
    "   b^{\\prime}(u^{-})))\\delta u = -A(u^{-})u^{-} + b(u^{-})$\n",
    "   with respect to $\\delta u$\n",
    "\n",
    "2. $u = u^{-} + \\omega\\delta u$\n",
    "\n",
    "3. $u^{-}\\ \\leftarrow\\ u$\n",
    "\n",
    "$\\gamma =1$ gives a Newton method while $\\gamma =0$ corresponds to\n",
    "Picard iteration.\n",
    "\n",
    "\n",
    "\n",
    "## Stopping criteria\n",
    "<div id=\"nonlin:systems:alg:terminate\"></div>\n",
    "\n",
    "\n",
    "Let $||\\cdot||$ be the standard Euclidean vector norm. Four termination\n",
    "criteria are much in use:\n",
    "\n",
    " * Absolute change in solution: $||u - u^{-}||\\leq \\epsilon_u$\n",
    "\n",
    " * Relative change in solution: $||u - u^{-}||\\leq \\epsilon_u ||u_0||$,\n",
    "   where $u_0$ denotes the start value of $u^{-}$ in the iteration\n",
    "\n",
    " * Absolute residual: $||F(u)|| \\leq \\epsilon_r$\n",
    "\n",
    " * Relative residual: $||F(u)|| \\leq \\epsilon_r ||F(u_0)||$\n",
    "\n",
    "To prevent divergent iterations to run forever,\n",
    "one terminates the iterations when\n",
    "the current number of iterations $k$ exceeds a maximum value $k_{\\max}$.\n",
    "\n",
    "For stationary problems, the relative criteria are most used since they are not sensitive to\n",
    "the characteristic size of $u$, which may depend on the underlying mesh and its resolution.  Nevertheless, the relative criteria\n",
    "can be misleading when the initial start value for the iteration is\n",
    "very close to the solution, since an unnecessary reduction in\n",
    "the error measure is enforced. \n",
    "For time-dependent problems, if the time-step is small then the previous solution may \n",
    "be a quite good guess for the unknown and in such cases the absolute criteria\n",
    "works better. It is common to combine the absolute and relative measures\n",
    "of the size of the residual,\n",
    "as in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto7\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "||F(u)|| \\leq \\epsilon_{rr} ||F(u_0)|| + \\epsilon_{ra},\n",
    "\\label{_auto7} \\tag{21}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\epsilon_{rr}$ is the tolerance in the relative criterion\n",
    "and $\\epsilon_{ra}$ is the tolerance in the absolute criterion.\n",
    "With a very good initial guess for the iteration\n",
    "(typically the solution of a differential\n",
    "equation at the previous time level), the term $||F(u_0)||$ is small\n",
    "and $\\epsilon_{ra}$ is the dominating tolerance. Otherwise,\n",
    "$\\epsilon_{rr} ||F(u_0)||$ and the relative criterion dominates.\n",
    "\n",
    "With the change in solution as criterion we can formulate a combined\n",
    "absolute and relative measure of the change in the solution:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto8\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "||\\delta u|| \\leq \\epsilon_{ur} ||u_0|| + \\epsilon_{ua},\n",
    "\\label{_auto8} \\tag{22}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ultimate termination criterion, combining the residual and\n",
    "the change in solution with a test on the maximum number\n",
    "of iterations, can be expressed as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto9\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "||F(u)|| \\leq \\epsilon_{rr} ||F(u_0)|| + \\epsilon_{ra}\n",
    "\\quad\\hbox{or}\\quad\n",
    "||\\delta u|| \\leq \\epsilon_{ur} ||u_0|| + \\epsilon_{ua}\n",
    "\\quad\\hbox{or}\\quad\n",
    "k>k_{\\max}{\\thinspace .}\n",
    "\\label{_auto9} \\tag{23}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonlinear PDEs\n",
    "\n",
    "# Linearization at the differential equation level\n",
    "<div id=\"nonlin:pdelevel\"></div>\n",
    "\n",
    "The attention is now turned to nonlinear partial differential\n",
    "equations (PDEs) and application of the techniques explained above for\n",
    "ODEs.  The model problem is a nonlinear diffusion equation for\n",
    "$u(\\boldsymbol{x},t)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:model:pde\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{\\partial u}{\\partial t} = \\nabla\\cdot ({\\alpha}(u)\\nabla u) + f(u),\\quad\n",
    "\\boldsymbol{x}\\in\\Omega,\\ t\\in (0,T],\n",
    "\\label{nonlin:pdelevel:model:pde} \\tag{1}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:model:Neumann\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "-{\\alpha}(u)\\frac{\\partial u}{\\partial n} = g,\\quad \\boldsymbol{x}\\in\\partial\\Omega_N,\\ \n",
    "t\\in (0,T],\n",
    "\\label{nonlin:pdelevel:model:Neumann} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:model:Dirichlet\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "u = u_0,\\quad \\boldsymbol{x}\\in\\partial\\Omega_D,\\ t\\in (0,T]{\\thinspace .}\n",
    "\\label{nonlin:pdelevel:model:Dirichlet} \\tag{3}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the present section, our aim is to discretize this problem in time\n",
    "and then present techniques for linearizing the time-discrete PDE\n",
    "problem \"at the PDE level\" such that we transform the nonlinear\n",
    "stationary PDE problem at each time level into a sequence of linear\n",
    "PDE problems, which can be solved using any method for linear\n",
    "PDEs. This strategy avoids the solution of systems of nonlinear\n",
    "algebraic equations.  In the section [1D stationary nonlinear differential equations](#nonlin:alglevel:1D) we shall take\n",
    "the opposite (and more common) approach: discretize the nonlinear\n",
    "problem in time and space first, and then solve the resulting\n",
    "nonlinear algebraic equations at each time level by the methods of\n",
    "the section \"Systems of nonlinear algebraic equations\" (for nonlinear ODEs). Very often, the two approaches are\n",
    "mathematically identical, so there is no preference from a\n",
    "computational efficiency point of view.  The details of the ideas\n",
    "sketched above will hopefully become clear through the forthcoming\n",
    "examples.\n",
    "\n",
    "\n",
    "## Explicit time integration\n",
    "<div id=\"nonlin:pdelevel:explicit\"></div>\n",
    "\n",
    "The nonlinearities in the PDE are trivial to deal with if we choose an\n",
    "explicit time integration method for ([1](#nonlin:pdelevel:model:pde)),\n",
    "such as the Forward Euler method:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_t^+ u = \\nabla\\cdot ({\\alpha}(u)\\nabla u) + f(u)]^n,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or written out,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n+1} - u^n}{\\Delta t} = \\nabla\\cdot ({\\alpha}(u^n)\\nabla u^n)\n",
    "+ f(u^n),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is a linear equation in the unknown $u^{n+1}$ with solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+1} = u^n + \\Delta t\\nabla\\cdot ({\\alpha}(u^n)\\nabla u^n) +\n",
    "\\Delta t f(u^n){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The disadvantage with this discretization is\n",
    "the strict stability criterion, e.g.,  $\\Delta t \\leq h^2/(6\\max\\alpha)$\n",
    "for the case $f=0$ and a standard 2nd-order finite difference discretization\n",
    "in 3D space with mesh cell sizes $h=\\Delta x=\\Delta y=\\Delta z$.\n",
    "\n",
    "<!-- BC -->\n",
    "\n",
    "## Backward Euler scheme and Picard iteration\n",
    "<div id=\"nonlin:pdelevel:Picard\"></div>\n",
    "\n",
    "A Backward Euler scheme for ([1](#nonlin:pdelevel:model:pde))\n",
    "reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_t^- u = \\nabla\\cdot ({\\alpha}(u)\\nabla u) + f(u)]^n{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Written out,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n} - u^{n-1}}{\\Delta t} = \\nabla\\cdot ({\\alpha}(u^n)\\nabla u^n)\n",
    "+ f(u^n){\\thinspace .}\n",
    "\\label{nonlin:pdelevel:pde:BE} \\tag{4}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a nonlinear PDE for the unknown function $u^n(\\boldsymbol{x})$. Such a\n",
    "PDE can be viewed as a time-independent PDE where\n",
    "$u^{n-1}(\\boldsymbol{x})$ is a known function.\n",
    "\n",
    "We introduce a Picard iteration with $k$ as iteration counter.\n",
    "A typical linearization of the $\\nabla\\cdot({\\alpha}(u^n)\\nabla u^n)$ term\n",
    "in iteration $k+1$ is to use the previously computed $u^{n,k}$\n",
    "approximation in the diffusion coefficient: ${\\alpha}(u^{n,k})$.\n",
    "The nonlinear source term is treated similarly: $f(u^{n,k})$.\n",
    "The unknown function $u^{n,k+1}$ then fulfills the linear PDE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Picard:k\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n,k+1} - u^{n-1}}{\\Delta t} = \\nabla\\cdot ({\\alpha}(u^{n,k})\n",
    "\\nabla u^{n,k+1})\n",
    "+ f(u^{n,k}){\\thinspace .}\n",
    "\\label{nonlin:pdelevel:pde:BE:Picard:k} \\tag{5}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial guess for the Picard iteration at this time level can be\n",
    "taken as the solution at the previous time level: $u^{n,0}=u^{n-1}$.\n",
    "\n",
    "We can alternatively apply the implementation-friendly\n",
    "notation where $u$ corresponds to\n",
    "the unknown we want to solve for, i.e., $u^{n,k+1}$ above, and $u^{-}$\n",
    "is the most recently computed value, $u^{n,k}$ above. Moreover,\n",
    "$u^{(1)}$ denotes the unknown function at the previous time level, $u^{n-1}$\n",
    "above. The PDE to be solved in a Picard iteration then looks like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Picard\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u - u^{(1)}}{\\Delta t} = \\nabla\\cdot ({\\alpha}(u^{-})\n",
    "\\nabla u)\n",
    "+ f(u^{-}){\\thinspace .}\n",
    "\\label{nonlin:pdelevel:pde:BE:Picard} \\tag{6}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the beginning of the iteration we start with the value from the\n",
    "previous time level: $u^{-}=u^{(1)}$, and after each\n",
    "iteration, $u^{-}$ is updated to $u$.\n",
    "\n",
    "**Remark on notation.**\n",
    "\n",
    "The previous derivations of the numerical scheme for time discretizations\n",
    "of PDEs have, strictly\n",
    "speaking, a somewhat sloppy notation, but it is much used and convenient\n",
    "to read. A more precise notation must\n",
    "distinguish clearly between the exact solution of the PDE problem,\n",
    "here denoted ${u_{\\small\\mbox{e}}}(\\boldsymbol{x},t)$, and the exact solution of the spatial\n",
    "problem, arising after time discretization at each time level,\n",
    "where ([4](#nonlin:pdelevel:pde:BE)) is an example. The latter\n",
    "is here represented as $u^n(\\boldsymbol{x})$ and is an approximation to\n",
    "${u_{\\small\\mbox{e}}}(\\boldsymbol{x},t_n)$. Then we have another approximation $u^{n,k}(\\boldsymbol{x})$\n",
    "to $u^n(\\boldsymbol{x})$ when solving the nonlinear PDE problem for\n",
    "$u^n$ by iteration methods, as in ([5](#nonlin:pdelevel:pde:BE:Picard:k)).\n",
    "\n",
    "In our notation, $u$ is a synonym for $u^{n,k+1}$ and $u^{(1)}$ is\n",
    "a synonym for $u^{n-1}$, inspired by what are natural variable names\n",
    "in a code.\n",
    "We will usually state the PDE problem in terms of $u$ and\n",
    "quickly redefine the symbol $u$ to mean the numerical approximation,\n",
    "while ${u_{\\small\\mbox{e}}}$ is not explicitly introduced unless we need to talk about\n",
    "the exact solution and the approximate solution at the same time.\n",
    "\n",
    "\n",
    "\n",
    "## Backward Euler scheme and Newton's method\n",
    "<div id=\"nonlin:pdelevel:Newton\"></div>\n",
    "\n",
    "\n",
    "At time level $n$, we have to solve the stationary PDE\n",
    "([4](#nonlin:pdelevel:pde:BE)). In the previous section, we\n",
    "saw how this can be done with Picard iterations.\n",
    "Another alternative is to apply the idea of Newton's method\n",
    "in a clever way.\n",
    "Normally, Newton's method is defined for systems of *algebraic equations*,\n",
    "but the idea of the method can be applied at the PDE level too.\n",
    "\n",
    "### Linearization via Taylor expansions\n",
    "\n",
    "Let $u^{n,k}$ be an approximation to the unknown $u^n$. We seek a\n",
    "better approximation on\n",
    "the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:Newton:ansatz\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u^{n} = u^{n,k} + \\delta u{\\thinspace .}\n",
    "\\label{nonlin:pdelevel:Newton:ansatz} \\tag{7}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea is to insert ([7](#nonlin:pdelevel:Newton:ansatz)) in\n",
    "([4](#nonlin:pdelevel:pde:BE)), Taylor expand the nonlinearities\n",
    "and keep only the terms that are\n",
    "linear in $\\delta u$ (which makes ([7](#nonlin:pdelevel:Newton:ansatz))\n",
    "an approximation for $u^{n}$). Then we can solve a linear PDE for\n",
    "the correction $\\delta u$ and use ([7](#nonlin:pdelevel:Newton:ansatz))\n",
    "to find a new approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n,k+1}=u^{n,k}+\\delta u\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to $u^{n}$.\n",
    "Repeating this procedure gives a sequence $u^{n,k+1}$, $k=0,1,\\ldots$\n",
    "that hopefully converges to the goal $u^n$.\n",
    "\n",
    "Let us carry out all the mathematical details for the nonlinear diffusion\n",
    "PDE discretized by the Backward Euler method.\n",
    "Inserting ([7](#nonlin:pdelevel:Newton:ansatz)) in\n",
    "([4](#nonlin:pdelevel:pde:BE)) gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Newton1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n,k} +\\delta u - u^{n-1}}{\\Delta t} =\n",
    "\\nabla\\cdot ({\\alpha}(u^{n,k} + \\delta u)\\nabla (u^{n,k}+\\delta u))\n",
    "+ f(u^{n,k}+\\delta u){\\thinspace .}\n",
    "\\label{nonlin:pdelevel:pde:BE:Newton1} \\tag{8}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can Taylor expand ${\\alpha}(u^{n,k} + \\delta u)$ and\n",
    "$f(u^{n,k}+\\delta u)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "{\\alpha}(u^{n,k} + \\delta u) & = {\\alpha}(u^{n,k}) + \\frac{d{\\alpha}}{du}(u^{n,k})\n",
    "\\delta u + {\\mathcal{O}(\\delta u^2)}\\approx {\\alpha}(u^{n,k}) + {\\alpha}^{\\prime}(u^{n,k})\\delta u,\\\\ \n",
    "f(u^{n,k}+\\delta u) &=  f(u^{n,k}) + \\frac{df}{du}(u^{n,k})\\delta u\n",
    "+ {\\mathcal{O}(\\delta u^2)}\\approx f(u^{n,k}) + f^{\\prime}(u^{n,k})\\delta u{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inserting the linear approximations of ${\\alpha}$ and $f$ in\n",
    "([8](#nonlin:pdelevel:pde:BE:Newton1)) results in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n,k} +\\delta u - u^{n-1}}{\\Delta t} =\n",
    "\\nabla\\cdot ({\\alpha}(u^{n,k})\\nabla u^{n,k}) + f(u^{n,k}) + \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\qquad \\nabla\\cdot ({\\alpha}(u^{n,k})\\nabla \\delta u)\n",
    "+ \\nabla\\cdot ({\\alpha}^{\\prime}(u^{n,k})\\delta u\\nabla u^{n,k}) + \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Newton2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\qquad \\nabla\\cdot ({\\alpha}^{\\prime}(u^{n,k})\\delta u\\nabla \\delta u)\n",
    "+ f^{\\prime}(u^{n,k})\\delta u{\\thinspace .}\n",
    "\\label{nonlin:pdelevel:pde:BE:Newton2} \\tag{9}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The term ${\\alpha}^{\\prime}(u^{n,k})\\delta u\\nabla \\delta u$ is of\n",
    "order $\\delta u^2$\n",
    "and is therefore omitted since we expect the correction $\\delta u$\n",
    "to be small ($\\delta u \\gg \\delta u^2$).\n",
    "Reorganizing the equation gives a PDE\n",
    "for $\\delta u$ that we can write in short form as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\delta F(\\delta u; u^{n,k}) = -F(u^{n,k}),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Newton2:F\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F(u^{n,k}) = \\frac{u^{n,k} - u^{n-1}}{\\Delta t} -\n",
    "\\nabla\\cdot ({\\alpha}(u^{n,k})\\nabla u^{n,k}) + f(u^{n,k}),\n",
    "\\label{nonlin:pdelevel:pde:BE:Newton2:F} \\tag{10}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\delta F(\\delta u; u^{n,k}) =\n",
    "- \\frac{1}{\\Delta t}\\delta u +\n",
    "\\nabla\\cdot ({\\alpha}(u^{n,k})\\nabla \\delta u) + \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\quad \\nabla\\cdot ({\\alpha}^{\\prime}(u^{n,k})\\delta u\\nabla u^{n,k})\n",
    "+ f^{\\prime}(u^{n,k})\\delta u{\\thinspace .}\n",
    "\\label{_auto1} \\tag{11}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $\\delta F$ is a linear function of $\\delta u$, and\n",
    "$F$ contains only terms that are known, such that\n",
    "the PDE for $\\delta u$ is indeed linear.\n",
    "\n",
    "**Observations.**\n",
    "\n",
    "The notational form $\\delta F = -F$ resembles the Newton system $J\\delta u =-F$\n",
    "for systems of algebraic equations, with $\\delta F$ as $J\\delta u$.\n",
    "The unknown vector in a linear system of algebraic equations enters\n",
    "the system as a linear operator in terms of a\n",
    "matrix-vector product ($J\\delta u$), while at\n",
    "the PDE level we have a linear differential operator instead\n",
    "($\\delta F$).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Similarity with Picard iteration\n",
    "\n",
    "We can rewrite the PDE for $\\delta u$ in a slightly different way too\n",
    "if we define $u^{n,k} + \\delta u$ as $u^{n,k+1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{n,k+1} - u^{n-1}}{\\Delta t} =\n",
    "\\nabla\\cdot ({\\alpha}(u^{n,k})\\nabla u^{n,k+1}) + f(u^{n,k})\\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\qquad  + \\nabla\\cdot ({\\alpha}^{\\prime}(u^{n,k})\\delta u\\nabla u^{n,k})\n",
    "+ f^{\\prime}(u^{n,k})\\delta u{\\thinspace .}\n",
    "\\label{_auto2} \\tag{12}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the first line is the same PDE as arises in the Picard\n",
    "iteration, while the remaining terms arise from the differentiations\n",
    "that are an inherent ingredient in Newton's method.\n",
    "\n",
    "### Implementation\n",
    "\n",
    "For coding we want to introduce $u$ for $u^n$, $u^{-}$ for $u^{n,k}$ and\n",
    "$u^{(1)}$ for $u^{n-1}$. The formulas for $F$ and $\\delta F$\n",
    "are then more clearly written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:pdelevel:pde:BE:Newton2:F2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F(u^{-}) = \\frac{u^{-} - u^{(1)}}{\\Delta t} -\n",
    "\\nabla\\cdot ({\\alpha}(u^{-})\\nabla u^{-}) + f(u^{-}),\n",
    "\\label{nonlin:pdelevel:pde:BE:Newton2:F2} \\tag{13}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\delta F(\\delta u; u^{-}) =\n",
    "- \\frac{1}{\\Delta t}\\delta u +\n",
    "\\nabla\\cdot ({\\alpha}(u^{-})\\nabla \\delta u) + \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\quad \\nabla\\cdot ({\\alpha}^{\\prime}(u^{-})\\delta u\\nabla u^{-})\n",
    "+ f^{\\prime}(u^{-})\\delta u{\\thinspace .}\n",
    "\\label{_auto3} \\tag{14}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The form that orders the PDE as the Picard iteration terms plus\n",
    "the Newton method's derivative terms becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u - u^{(1)}}{\\Delta t} =\n",
    "\\nabla\\cdot ({\\alpha}(u^{-})\\nabla u) + f(u^{-}) + \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\qquad  \\gamma(\\nabla\\cdot ({\\alpha}^{\\prime}(u^{-})(u - u^{-})\\nabla u^{-})\n",
    "+ f^{\\prime}(u^{-})(u - u^{-})){\\thinspace .}\n",
    "\\label{_auto4} \\tag{15}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Picard and full Newton versions correspond to\n",
    "$\\gamma=0$ and $\\gamma=1$, respectively.\n",
    "\n",
    "\n",
    "\n",
    "### Derivation with alternative notation\n",
    "\n",
    "Some may prefer to derive the linearized PDE for $\\delta u$ using\n",
    "the more compact notation. We start with inserting $u^n=u^{-}+\\delta u$\n",
    "to get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{u^{-} +\\delta u - u^{n-1}}{\\Delta t} =\n",
    "\\nabla\\cdot ({\\alpha}(u^{-} + \\delta u)\\nabla (u^{-}+\\delta u))\n",
    "+ f(u^{-}+\\delta u){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taylor expanding,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "{\\alpha}(u^{-} + \\delta u) & \\approx {\\alpha}(u^{-}) + {\\alpha}^{\\prime}(u^{-})\\delta u,\\\\ \n",
    "f(u^{-}+\\delta u) & \\approx f(u^{-}) + f^{\\prime}(u^{-})\\delta u,\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and inserting these expressions gives a less cluttered PDE for $\\delta u$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{u^{-} +\\delta u - u^{n-1}}{\\Delta t} &=\n",
    "\\nabla\\cdot ({\\alpha}(u^{-})\\nabla u^{-}) + f(u^{-}) + \\\\ \n",
    "&\\qquad \\nabla\\cdot ({\\alpha}(u^{-})\\nabla \\delta u)\n",
    "+ \\nabla\\cdot ({\\alpha}^{\\prime}(u^{-})\\delta u\\nabla u^{-}) + \\\\ \n",
    "&\\qquad \\nabla\\cdot ({\\alpha}^{\\prime}(u^{-})\\delta u\\nabla \\delta u)\n",
    "+ f^{\\prime}(u^{-})\\delta u{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crank-Nicolson discretization\n",
    "<div id=\"nonlin:pdelevel:Picard:CN\"></div>\n",
    "\n",
    "A Crank-Nicolson discretization of\n",
    "([1](#nonlin:pdelevel:model:pde)) applies a centered difference\n",
    "at $t_{n+\\frac{1}{2}}$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_t u = \\nabla\\cdot ({\\alpha}(u)\\nabla u) + f(u)]^{n+\\frac{1}{2}}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard technique is to apply an arithmetic average for\n",
    "quantities defined between two mesh points, e.g.,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+\\frac{1}{2}}\\approx \\frac{1}{2}(u^n + u^{n+1}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, with nonlinear terms we have many choices of formulating\n",
    "an arithmetic mean:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto5\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "[f(u)]^{n+\\frac{1}{2}} \\approx f(\\frac{1}{2}(u^n + u^{n+1}))\n",
    "= [f(\\overline{u}^t)]^{n+\\frac{1}{2}},\n",
    "\\label{_auto5} \\tag{16}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto6\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "[f(u)]^{n+\\frac{1}{2}} \\approx \\frac{1}{2}(f(u^n) + f(u^{n+1}))\n",
    "=[\\overline{f(u)}^t]^{n+\\frac{1}{2}},\n",
    "\\label{_auto6} \\tag{17}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto7\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "[{\\alpha}(u)\\nabla u]^{n+\\frac{1}{2}} \\approx\n",
    "{\\alpha}(\\frac{1}{2}(u^n + u^{n+1}))\\nabla (\\frac{1}{2}(u^n + u^{n+1}))\n",
    "= [{\\alpha}(\\overline{u}^t)\\nabla \\overline{u}^t]^{n+\\frac{1}{2}},\n",
    "\\label{_auto7} \\tag{18}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto8\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "[{\\alpha}(u)\\nabla u]^{n+\\frac{1}{2}} \\approx\n",
    "\\frac{1}{2}({\\alpha}(u^n) + {\\alpha}(u^{n+1}))\\nabla (\\frac{1}{2}(u^n + u^{n+1}))\n",
    "= [\\overline{{\\alpha}(u)}^t\\nabla\\overline{u}^t]^{n+\\frac{1}{2}},\n",
    "\\label{_auto8} \\tag{19}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto9\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "[{\\alpha}(u)\\nabla u]^{n+\\frac{1}{2}} \\approx\n",
    "\\frac{1}{2}({\\alpha}(u^n)\\nabla u^n + {\\alpha}(u^{n+1})\\nabla u^{n+1})\n",
    "= [\\overline{{\\alpha}(u)\\nabla u}^t]^{n+\\frac{1}{2}}{\\thinspace .}\n",
    "\\label{_auto9} \\tag{20}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A big question is whether there are significant differences in accuracy\n",
    "between taking the products of arithmetic means or taking the arithmetic\n",
    "mean of products. [nonlin:exer:products:arith:mean](#nonlin:exer:products:arith:mean) investigates\n",
    "this question, and the answer is that the approximation is\n",
    "${\\mathcal{O}(\\Delta t^2)}$ in both cases. \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# 1D stationary nonlinear differential equations\n",
    "<div id=\"nonlin:alglevel:1D\"></div>\n",
    "\n",
    "The section [Linearization at the differential equation level](#nonlin:pdelevel) presented methods for linearizing\n",
    "time-discrete PDEs directly prior to discretization in space.  We can\n",
    "alternatively carry out the discretization in space of the\n",
    "time-discrete nonlinear PDE problem and get a system of nonlinear\n",
    "algebraic equations, which can be solved by Picard iteration or\n",
    "Newton's method as presented in the section \"Systems of nonlinear algebraic equations\" (for nonlinear ODEs).\n",
    "This latter approach will now be described in detail.\n",
    "\n",
    "We shall work with the 1D problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:pde\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-({\\alpha}(u)u^{\\prime})^{\\prime} + au = f(u),\\quad x\\in (0,L),\n",
    "\\quad {\\alpha}(u(0))u^{\\prime}(0) = C,\\ u(L)=D\n",
    "{\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:pde} \\tag{21}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem ([21](#nonlin:alglevel:1D:pde)) arises from the stationary\n",
    "limit of a diffusion equation,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:pde:tver\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{\\partial u}{\\partial t} = \\frac{\\partial}{\\partial x}\\left(\n",
    "\\alpha(u)\\frac{\\partial u}{\\partial x}\\right) - au + f(u),\n",
    "\\label{nonlin:alglevel:1D:pde:tver} \\tag{22}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as $t\\rightarrow\\infty$ and $\\partial u/\\partial t\\rightarrow 0$.\n",
    "Alternatively, the problem ([21](#nonlin:alglevel:1D:pde)) arises\n",
    "at each time level from implicit time discretization of\n",
    "([22](#nonlin:alglevel:1D:pde:tver)). For example, a Backward Euler\n",
    "scheme for ([22](#nonlin:alglevel:1D:pde:tver)) leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:pde:tver:BE\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{u^{n}-u^{n-1}}{\\Delta t} =\n",
    "\\frac{d}{dx}\\left(\n",
    "\\alpha(u^n)\\frac{du^n}{dx}\\right) - au^n + f(u^n){\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:pde:tver:BE} \\tag{23}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Introducing $u(x)$ for $u^n(x)$, $u^{(1)}$ for $u^{n-1}$, and defining $f(u)$\n",
    "in ([21](#nonlin:alglevel:1D:pde)) to be $f(u)$ in\n",
    "([23](#nonlin:alglevel:1D:pde:tver:BE)) plus $u^{n-1}/\\Delta t$, gives\n",
    "([21](#nonlin:alglevel:1D:pde)) with $a=1/\\Delta t$.\n",
    "\n",
    "\n",
    "## Finite difference discretization\n",
    "<div id=\"nonlin:alglevel:1D:fd\"></div>\n",
    "\n",
    "\n",
    "The nonlinearity in the differential equation\n",
    "([21](#nonlin:alglevel:1D:pde)) poses no more difficulty than a variable\n",
    "coefficient, as in the term $({\\alpha}(x)u^{\\prime})^{\\prime}$.  We can\n",
    "therefore use a standard finite difference approach to discretizing\n",
    "the Laplace term with a variable coefficient:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[-D_x{\\alpha} D_x u +au = f]_i{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Writing this out for a uniform mesh with points $x_i=i\\Delta x$,\n",
    "$i=0,\\ldots,N_x$, leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:deq0\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{1}{\\Delta x^2}\n",
    "\\left({\\alpha}_{i+\\frac{1}{2}}(u_{i+1}-u_i) -\n",
    "{\\alpha}_{i-\\frac{1}{2}}(u_{i}-u_{i-1})\\right)\n",
    "+ au_i = f(u_i){\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:fd:deq0} \\tag{24}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This equation is valid at all the mesh points $i=0,1,\\ldots,N_x-1$.\n",
    "At $i=N_x$ we have the Dirichlet condition $u_i=D$.\n",
    "The only difference from the case with $({\\alpha}(x)u^{\\prime})^{\\prime}$ and $f(x)$ is that\n",
    "now ${\\alpha}$ and $f$ are functions of $u$ and not only of $x$:\n",
    "$({\\alpha}(u(x))u^{\\prime})^{\\prime}$ and $f(u(x))$.\n",
    "\n",
    "The quantity ${\\alpha}_{i+\\frac{1}{2}}$, evaluated between two mesh points,\n",
    "needs a comment. Since ${\\alpha}$ depends on $u$ and $u$ is only known\n",
    "at the mesh points, we need to express ${\\alpha}_{i+\\frac{1}{2}}$ in\n",
    "terms of $u_i$ and $u_{i+1}$. For this purpose we use an arithmetic\n",
    "mean, although a harmonic mean is also common in this context if\n",
    "${\\alpha}$ features large jumps.\n",
    "There are two choices of arithmetic means:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:dfc:mean:u\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\alpha}_{i+\\frac{1}{2}} \\approx\n",
    "{\\alpha}(\\frac{1}{2}(u_i + u_{i+1}) =\n",
    "[{\\alpha}(\\overline{u}^x)]^{i+\\frac{1}{2}},\n",
    "\\label{nonlin:alglevel:1D:fd:dfc:mean:u} \\tag{25}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:dfc:mean:dfc\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "{\\alpha}_{i+\\frac{1}{2}} \\approx\n",
    "\\frac{1}{2}({\\alpha}(u_i) + {\\alpha}(u_{i+1})) = [\\overline{{\\alpha}(u)}^x]^{i+\\frac{1}{2}}\n",
    "\\label{nonlin:alglevel:1D:fd:dfc:mean:dfc} \\tag{26}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation ([24](#nonlin:alglevel:1D:fd:deq0)) with\n",
    "the latter approximation then looks like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\frac{1}{2\\Delta x^2}\n",
    "\\left(({\\alpha}(u_i)+{\\alpha}(u_{i+1}))(u_{i+1}-u_i) -\n",
    "({\\alpha}(u_{i-1})+{\\alpha}(u_{i}))(u_{i}-u_{i-1})\\right)\\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:deq\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\qquad\\qquad + au_i = f(u_i),\n",
    "\\label{nonlin:alglevel:1D:fd:deq} \\tag{27}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or written more compactly,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[-D_x\\overline{{\\alpha}}^x D_x u +au = f]_i{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At mesh point $i=0$ we have the boundary condition ${\\alpha}(u)u^{\\prime}=C$,\n",
    "which is discretized by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[{\\alpha}(u)D_{2x}u = C]_0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "meaning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:Neumann:x0\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\alpha}(u_0)\\frac{u_{1} - u_{-1}}{2\\Delta x} = C{\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:fd:Neumann:x0} \\tag{28}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fictitious value $u_{-1}$ can be eliminated with the aid\n",
    "of ([27](#nonlin:alglevel:1D:fd:deq)) for $i=0$.\n",
    "Formally, ([27](#nonlin:alglevel:1D:fd:deq)) should be solved with\n",
    "respect to $u_{i-1}$ and that value (for $i=0$) should be inserted in\n",
    "([28](#nonlin:alglevel:1D:fd:Neumann:x0)), but it is algebraically\n",
    "much easier to do it the other way around. Alternatively, one can\n",
    "use a ghost cell $[-\\Delta x,0]$ and update the $u_{-1}$ value\n",
    "in the ghost cell according to ([28](#nonlin:alglevel:1D:fd:Neumann:x0))\n",
    "after every Picard or Newton iteration. Such an approach means that\n",
    "we use a known $u_{-1}$ value in ([27](#nonlin:alglevel:1D:fd:deq))\n",
    "from the previous iteration.\n",
    "\n",
    "## Solution of algebraic equations\n",
    "\n",
    "### The structure of the equation system\n",
    "\n",
    "The nonlinear algebraic equations ([27](#nonlin:alglevel:1D:fd:deq)) are\n",
    "of the form $A(u)u = b(u)$ with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{i,i} &= \\frac{1}{2\\Delta x^2}({\\alpha}(u_{i-1}) + 2{\\alpha}(u_{i})\n",
    "+ {\\alpha}(u_{i+1})) + a,\\\\ \n",
    "A_{i,i-1} &= -\\frac{1}{2\\Delta x^2}({\\alpha}(u_{i-1}) + {\\alpha}(u_{i})),\\\\ \n",
    "A_{i,i+1} &= -\\frac{1}{2\\Delta x^2}({\\alpha}(u_{i}) + {\\alpha}(u_{i+1})),\\\\ \n",
    "b_i &= f(u_i){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix $A(u)$ is tridiagonal: $A_{i,j}=0$ for $j > i+1$ and $j < i-1$.\n",
    "\n",
    "The above expressions are valid for internal mesh points $1\\leq i\\leq N_x-1$.\n",
    "For $i=0$ we need to express $u_{i-1}=u_{-1}$ in terms of $u_1$ using\n",
    "([28](#nonlin:alglevel:1D:fd:Neumann:x0)):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:Neumann:x0:um1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u_{-1} = u_1 -\\frac{2\\Delta x}{{\\alpha}(u_0)}C{\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:fd:Neumann:x0:um1} \\tag{29}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This value must be inserted in $A_{0,0}$. The expression for $A_{i,i+1}$\n",
    "applies for $i=0$, and $A_{i,i-1}$ does not enter the system when $i=0$.\n",
    "\n",
    "Regarding the last equation, its form depends on whether we include\n",
    "the Dirichlet condition $u(L)=D$, meaning $u_{N_x}=D$, in the\n",
    "nonlinear algebraic equation system or not. Suppose we choose\n",
    "$(u_0,u_1,\\ldots,u_{N_x-1})$ as unknowns, later referred to as\n",
    "*systems without Dirichlet conditions*. The last equation\n",
    "corresponds to $i=N_x-1$. It involves the boundary value $u_{N_x}$,\n",
    "which is substituted by $D$. If the unknown vector includes the\n",
    "boundary value, $(u_0,u_1,\\ldots,u_{N_x})$, later referred to as\n",
    "*system including Dirichlet conditions*, the equation for $i=N_x-1$\n",
    "just involves the unknown $u_{N_x}$, and the final equation becomes\n",
    "$u_{N_x}=D$, corresponding to $A_{i,i}=1$ and $b_i=D$ for $i=N_x$.\n",
    "\n",
    "### Picard iteration\n",
    "\n",
    "The obvious Picard iteration scheme is to use previously computed\n",
    "values of $u_i$ in $A(u)$ and $b(u)$, as described more in detail in\n",
    "the section [nonlin:systems:alg](#nonlin:systems:alg). With the notation $u^{-}$ for the\n",
    "most recently computed value of $u$, we have the system $F(u)\\approx\n",
    "\\hat F(u) = A(u^{-})u - b(u^{-})$, with $F=(F_0,F_1,\\ldots,F_m)$,\n",
    "$u=(u_0,u_1,\\ldots,u_m)$.  The index $m$ is $N_x$ if the system\n",
    "includes the Dirichlet condition as a separate equation and $N_x-1$\n",
    "otherwise.  The matrix $A(u^{-})$ is tridiagonal, so the solution\n",
    "procedure is to fill a tridiagonal matrix data structure and the\n",
    "right-hand side vector with the right numbers and call a Gaussian\n",
    "elimination routine for tridiagonal linear systems.\n",
    "\n",
    "### Mesh with two cells\n",
    "\n",
    "It helps on the understanding of the details to write out all the\n",
    "mathematics in a specific\n",
    "case with a small mesh, say just two cells ($N_x=2$). We use $u^{-}_i$\n",
    "for the $i$-th component in $u^{-}$.\n",
    "\n",
    "The starting point is the basic expressions for the\n",
    "nonlinear equations at mesh point $i=0$ and $i=1$ are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:2x2:x0\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "A_{0,-1}u_{-1} + A_{0,0}u_0 + A_{0,1}u_1 = b_0,\n",
    "\\label{nonlin:alglevel:1D:fd:2x2:x0} \\tag{30}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:2x2:x1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "A_{1,0}u_{0} + A_{1,1}u_1 + A_{1,2}u_2 = b_1{\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:fd:2x2:x1} \\tag{31}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation ([30](#nonlin:alglevel:1D:fd:2x2:x0)) written out reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{1}{2\\Delta x^2}(& -({\\alpha}(u_{-1}) + {\\alpha}(u_{0}))u_{-1}\\, +\\\\ \n",
    "& ({\\alpha}(u_{-1}) + 2{\\alpha}(u_{0}) + {\\alpha}(u_{1}))u_0\\, -\\\\ \n",
    "& ({\\alpha}(u_{0}) + {\\alpha}(u_{1})))u_1 + au_0\n",
    "=f(u_0){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must then replace $u_{-1}$ by\n",
    "([29](#nonlin:alglevel:1D:fd:Neumann:x0:um1)).\n",
    "With Picard iteration we get\n",
    "<!-- u_{-1} = u_1 -\\frac{2\\Delta x}{{\\alpha}(u_0)}C -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{1}{2\\Delta x^2}(& -({\\alpha}(u^-_{-1}) + 2{\\alpha}(u^-_{0})\n",
    "+ {\\alpha}(u^-_{1}))u_1\\, +\\\\ \n",
    "&({\\alpha}(u^-_{-1}) + 2{\\alpha}(u^-_{0}) + {\\alpha}(u^-_{1}))u_0\n",
    " + au_0\\\\ \n",
    "&=f(u^-_0) -\n",
    "\\frac{1}{{\\alpha}(u^-_0)\\Delta x}({\\alpha}(u^-_{-1}) + {\\alpha}(u^-_{0}))C,\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^-_{-1} = u_1^- -\\frac{2\\Delta x}{{\\alpha}(u^-_0)}C{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation ([31](#nonlin:alglevel:1D:fd:2x2:x1)) contains the unknown $u_2$\n",
    "for which we have a Dirichlet condition. In case we omit the\n",
    "condition as a separate equation, ([31](#nonlin:alglevel:1D:fd:2x2:x1))\n",
    "with Picard iteration becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{1}{2\\Delta x^2}(&-({\\alpha}(u^-_{0}) + {\\alpha}(u^-_{1}))u_{0}\\, + \\\\ \n",
    "&({\\alpha}(u^-_{0}) + 2{\\alpha}(u^-_{1}) + {\\alpha}(u^-_{2}))u_1\\, -\\\\ \n",
    "&({\\alpha}(u^-_{1}) + {\\alpha}(u^-_{2})))u_2 + au_1\n",
    "=f(u^-_1){\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must now move the $u_2$ term to the right-hand side and replace all\n",
    "occurrences of $u_2$ by $D$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{1}{2\\Delta x^2}(&-({\\alpha}(u^-_{0}) + {\\alpha}(u^-_{1}))u_{0}\\, +\\\\ \n",
    "& ({\\alpha}(u^-_{0}) + 2{\\alpha}(u^-_{1}) + {\\alpha}(D))u_1 + au_1\\\\ \n",
    "&=f(u^-_1) + \\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{1}) + {\\alpha}(D))D{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two equations can be written as a $2\\times 2$ system:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{cc}\n",
    "B_{0,0}& B_{0,1}\\\\ \n",
    "B_{1,0} & B_{1,1}\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "u_0\\\\ \n",
    "u_1\n",
    "\\end{array}\\right)\n",
    "=\n",
    "\\left(\\begin{array}{c}\n",
    "d_0\\\\ \n",
    "d_1\n",
    "\\end{array}\\right),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto10\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "B_{0,0} =\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{-1}) + 2{\\alpha}(u^-_{0}) + {\\alpha}(u^-_{1}))\n",
    "+ a\n",
    "\\label{_auto10} \\tag{32}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto11\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "B_{0,1} =\n",
    "-\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{-1}) + 2{\\alpha}(u^-_{0})\n",
    "+ {\\alpha}(u^-_{1})),\n",
    "\\label{_auto11} \\tag{33}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto12\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "B_{1,0} =\n",
    "-\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{0}) + {\\alpha}(u^-_{1})),\n",
    "\\label{_auto12} \\tag{34}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto13\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "B_{1,1} =\n",
    "\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{0}) + 2{\\alpha}(u^-_{1}) + {\\alpha}(D)) + a,\n",
    "\\label{_auto13} \\tag{35}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto14\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "d_0 =\n",
    "f(u^-_0) -\n",
    "\\frac{1}{{\\alpha}(u^-_0)\\Delta x}({\\alpha}(u^-_{-1}) + {\\alpha}(u^-_{0}))C,\n",
    "\\label{_auto14} \\tag{36}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto15\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "d_1 = f(u^-_1) + \\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{1}) + {\\alpha}(D))D{\\thinspace .}\n",
    "\\label{_auto15} \\tag{37}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system with the Dirichlet condition becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{ccc}\n",
    "B_{0,0}& B_{0,1} & 0\\\\ \n",
    "B_{1,0} & B_{1,1} & B_{1,2}\\\\ \n",
    "0  & 0 & 1\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "u_0\\\\ \n",
    "u_1\\\\ \n",
    "u_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "\\left(\\begin{array}{c}\n",
    "d_0\\\\ \n",
    "d_1\\\\ \n",
    "D\n",
    "\\end{array}\\right),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto16\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "B_{1,1} =\n",
    "\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{0}) + 2{\\alpha}(u^-_{1}) + {\\alpha}(u_2)) + a,\n",
    "\\label{_auto16} \\tag{38}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto17\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "B_{1,2} = -\n",
    "\\frac{1}{2\\Delta x^2}({\\alpha}(u^-_{1}) + {\\alpha}(u_2))),\n",
    "\\label{_auto17} \\tag{39}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto18\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "d_1 = f(u^-_1){\\thinspace .}\n",
    "\\label{_auto18} \\tag{40}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other entries are as in the $2\\times 2$ system.\n",
    "\n",
    "\n",
    "### Newton's method\n",
    "\n",
    "The Jacobian must be derived in order to use Newton's method. Here it means\n",
    "that we need to differentiate $F(u)=A(u)u - b(u)$ with respect to\n",
    "the unknown parameters\n",
    "$u_0,u_1,\\ldots,u_m$ ($m=N_x$ or $m=N_x-1$, depending on whether the\n",
    "Dirichlet condition is included in the nonlinear system $F(u)=0$ or not).\n",
    "Nonlinear equation number $i$ of ([27](#nonlin:alglevel:1D:fd:deq))\n",
    "has the structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_i = A_{i,i-1}(u_{i-1},u_i)u_{i-1} +\n",
    "A_{i,i}(u_{i-1},u_i,u_{i+1})u_i +\n",
    "A_{i,i+1}(u_i, u_{i+1})u_{i+1} - b_i(u_i){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the Jacobian requires careful differentiation. For example,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{\\partial}{\\partial u_i}(A_{i,i}(u_{i-1},u_i,u_{i+1})u_i) &=\n",
    "\\frac{\\partial A_{i,i}}{\\partial u_i}u_i + A_{i,i}\n",
    "\\frac{\\partial u_i}{\\partial u_i}\\\\ \n",
    "&=\n",
    "\\frac{\\partial}{\\partial u_i}(\n",
    "\\frac{1}{2\\Delta x^2}({\\alpha}(u_{i-1}) + 2{\\alpha}(u_{i})\n",
    "+{\\alpha}(u_{i+1})) + a)u_i +\\\\ \n",
    "&\\quad\\frac{1}{2\\Delta x^2}({\\alpha}(u_{i-1}) + 2{\\alpha}(u_{i})\n",
    "+{\\alpha}(u_{i+1})) + a\\\\ \n",
    "&= \\frac{1}{2\\Delta x^2}(2{\\alpha}^\\prime (u_i)u_i\n",
    "+{\\alpha}(u_{i-1}) + 2{\\alpha}(u_{i})\n",
    "+{\\alpha}(u_{i+1})) + a{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete Jacobian becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "J_{i,i} &= \\frac{\\partial F_i}{\\partial u_i}\n",
    "= \\frac{\\partial A_{i,i-1}}{\\partial u_i}u_{i-1}\n",
    "+ \\frac{\\partial A_{i,i}}{\\partial u_i}u_i\n",
    "+ A_{i,i}\n",
    "+ \\frac{\\partial A_{i,i+1}}{\\partial u_i}u_{i+1}\n",
    "- \\frac{\\partial b_i}{\\partial u_{i}}\\\\ \n",
    "&=\n",
    "\\frac{1}{2\\Delta x^2}(\n",
    "-{\\alpha}^{\\prime}(u_i)u_{i-1}\n",
    "+2{\\alpha}^{\\prime}(u_i)u_{i}\n",
    "+{\\alpha}(u_{i-1}) + 2{\\alpha}(u_i) + {\\alpha}(u_{i+1})) +\\\\ \n",
    "&\\quad a\n",
    "-\\frac{1}{2\\Delta x^2}{\\alpha}^{\\prime}(u_{i})u_{i+1}\n",
    "- b^{\\prime}(u_i),\\\\ \n",
    "J_{i,i-1} &= \\frac{\\partial F_i}{\\partial u_{i-1}}\n",
    "= \\frac{\\partial A_{i,i-1}}{\\partial u_{i-1}}u_{i-1}\n",
    "+ A_{i-1,i}\n",
    "+ \\frac{\\partial A_{i,i}}{\\partial u_{i-1}}u_i\n",
    "- \\frac{\\partial b_i}{\\partial u_{i-1}}\\\\ \n",
    "&=\n",
    "\\frac{1}{2\\Delta x^2}(\n",
    "-{\\alpha}^{\\prime}(u_{i-1})u_{i-1} - ({\\alpha}(u_{i-1}) + {\\alpha}(u_i))\n",
    "+ {\\alpha}^{\\prime}(u_{i-1})u_i),\\\\ \n",
    "J_{i,i+1} &= \\frac{\\partial A_{i,i+1}}{\\partial u_{i-1}}u_{i+1}\n",
    "+ A_{i+1,i} +\n",
    "\\frac{\\partial A_{i,i}}{\\partial u_{i+1}}u_i\n",
    "- \\frac{\\partial b_i}{\\partial u_{i+1}}\\\\ \n",
    "&=\\frac{1}{2\\Delta x^2}(\n",
    "-{\\alpha}^{\\prime}(u_{i+1})u_{i+1} - ({\\alpha}(u_{i}) + {\\alpha}(u_{i+1}))\n",
    "+ {\\alpha}^{\\prime}(u_{i+1})u_i)\n",
    "{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The explicit expression for nonlinear equation number $i$,\n",
    "$F_i(u_0,u_1,\\ldots)$, arises from moving the $f(u_i)$ term in\n",
    "([27](#nonlin:alglevel:1D:fd:deq)) to the left-hand side:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_i = -\\frac{1}{2\\Delta x^2}\n",
    "\\left(({\\alpha}(u_i)+{\\alpha}(u_{i+1}))(u_{i+1}-u_i) -\n",
    "({\\alpha}(u_{i-1})+{\\alpha}(u_{i}))(u_{i}-u_{i-1})\\right)\\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"nonlin:alglevel:1D:fd:deq2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}  \n",
    "\\qquad\\qquad + au_i - f(u_i) = 0{\\thinspace .}\n",
    "\\label{nonlin:alglevel:1D:fd:deq2} \\tag{41}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the boundary point $i=0$, $u_{-1}$ must be replaced using\n",
    "the formula ([29](#nonlin:alglevel:1D:fd:Neumann:x0:um1)).\n",
    "When the Dirichlet condition at $i=N_x$ is not a part of the\n",
    "equation system, the last equation $F_m=0$ for $m=N_x-1$\n",
    "involves the quantity $u_{N_x-1}$ which must be replaced by $D$.\n",
    "If $u_{N_x}$ is treated as an unknown in the system, the\n",
    "last equation $F_m=0$ has $m=N_x$ and reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{N_x}(u_0,\\ldots,u_{N_x}) = u_{N_x} - D = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar replacement of $u_{-1}$ and $u_{N_x}$ must be done in\n",
    "the Jacobian for the first and last row. When $u_{N_x}$\n",
    "is included as an unknown, the last row in the Jacobian\n",
    "must help implement the condition $\\delta u_{N_x}=0$, since\n",
    "we assume that $u$ contains the right Dirichlet value\n",
    "at the beginning of the iteration ($u_{N_x}=D$), and then\n",
    "the Newton update should be zero for $i=0$, i.e., $\\delta u_{N_x}=0$.\n",
    "This also forces the right-hand side to be $b_i=0$, $i=N_x$.\n",
    "\n",
    "We have seen, and can see from the present example, that the\n",
    "linear system in Newton's method contains all the terms present\n",
    "in the system that arises in the Picard iteration method.\n",
    "The extra terms in Newton's method can be multiplied by a factor\n",
    "such that it is easy to program one linear system and set this\n",
    "factor to 0 or 1 to generate the Picard or Newton system.\n",
    "\n",
    "<!-- Remark: Neumann cond at x=L and Dirichlet at x=0 leads to different -->\n",
    "<!-- numbering of unknowns and u at mesh points. Must address this -->\n",
    "<!-- in a remark and treat it properly in diffu. -->"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}