The attention is now turned to nonlinear partial differential equations (PDEs) and application of the techniques explained above for ODEs. The model problem is a nonlinear diffusion equation $$ \begin{align} \frac{\partial u}{\partial t} &= \nabla\cdot (\dfc(u)\nabla u) + f(u),\quad &\x\in\Omega,\ t\in (0,T], \tag{17} \\ -\dfc(u)\frac{\partial u}{\partial n} &= g,\quad &\x\in\partial\Omega_N,\ t\in (0,T], \tag{18} \\ u &= u_0,\quad &\x\in\partial\Omega_D,\ t\in (0,T]\tp \tag{19} \end{align} $$
Our aim is to discretize the problem in time and then present techniques for linearizing the time-discrete PDE problem "at the PDE level" such that we transform the nonlinear stationary PDE problems at each time level into a sequence of linear PDE problems, which can be solved using any method for linear PDEs. This strategy avoids the solution systems of nonlinear algebraic equations. In the section Discretization of 1D stationary nonlinear differential equations we shall take the opposite (and more common) approach: discretize the nonlinear problem in time and space first, and then solve the resulting nonlinear algebraic equations at each time level by the methods of the section Systems of nonlinear algebraic equations.
The nonlinearities in the PDE are trivial to deal with if we choose an explicit time integration method for (17), such as the Forward Euler method: $$ [D_t^+ u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^n,$$ or written out, $$ \frac{u^{n+1} - u^n}{\Delta t} = \nabla\cdot (\dfc(u^n)\nabla u^n) + f(u^n),$$ which is a linear equation in the unknown \( u^{n+1} \) with solution $$ u^{n+1} = u^n + \Delta t\nabla\cdot (\dfc(u^n)\nabla u^n) + \Delta t f(u^n)\tp $$
The disadvantage with this discretization is usually thought to be the stability criterion $$ \Delta t \leq \frac{1}{\max\alpha}(\Delta x^2 + \Delta y^2 + \Delta z^2),$$ for the case \( f=0 \) and a standard 2nd-order finite difference discretization in space with mesh cell sizes \( \Delta x \), \( \Delta y \), and \( \Delta z \) in the various spatial directions.
A Backward Euler scheme for (17) reads $$ [D_t^- u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^n\tp$$ Written out, $$ \begin{equation} \frac{u^{n} - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^n)\nabla u^n) + f(u^n)\tp \tag{20} \end{equation} $$ This is a nonlinear PDE for the unknown function \( u^n(\x) \). Such a PDE can be viewed as a time-independent PDE where \( u^{n-1}(\x) \) is a known function.
We introduce a Picard iteration with \( k \) as iteration counter. A typical linearization of the \( \nabla\cdot\dfc(u^n)\nabla u^n \) term in iteration \( k+1 \) is to use the previously computed \( u^{n,k} \) approximation in the diffusion coefficient: \( \dfc(u^{n,k}) \). The nonlinear source term is treated similarly: \( f(u^{n,k}) \). The unknown function \( u^{n,k+1} \) then fulfills the linear PDE $$ \begin{equation} \frac{u^{n,k+1} - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^{n,k}) \nabla u^{n,k+1}) + f(u^{n,k})\tp \tag{21} \end{equation} $$ The initial guess for the Picard iteration at this time level can be taken as the solution at the previous time level: \( u^{n,0}=u^{n-1} \).
We can alternatively apply the implementation-friendly notation where \( u \) corresponds to the unknown we want to solve for, i.e., \( u^{n,k+1} \) above, and \( u^{-} \) is the most recently computed value, \( u^{n,k} \) above. Moreover, \( u^{(1)} \) denotes the unknown function at the previous time level, \( u^{n-1} \) above. The PDE to be solved in a Picard iteration then looks like $$ \begin{equation} \frac{u - u^{(1)}}{\Delta t} = \nabla\cdot (\dfc(u^{-}) \nabla u) + f(u^{-})\tp \tag{22} \end{equation} $$ At the beginning of the iteration we start with the value from the previous time level: \( u^{-}=u^{(1)} \), and after each iteration, \( u^{-} \) is updated to \( u \).
However, we will usually state the PDE problem in terms of \( u \) and quickly redefine the symbol \( u \) to mean the numerical approximation, while \( \uex \) is not explicitly introduced unless we need to talk about the exact solution and the approximate solution at the same time.
At time level \( n \) we have to solve the stationary PDE (20), this time with Newton's method. Normally, Newton's method is defined for systems of algebraic equations, but the idea of the method can be applied at the PDE level too.
Let \( u^{n,k} \) be an approximation to the unknown \( u^n \). We seek a better approximation on the form $$ \begin{equation} u^{n} = u^{n,k} + \delta u\tp \tag{23} \end{equation} $$ The idea is to insert (23) in (20), Taylor expand the nonlinearities and keep only the terms that are linear in \( \delta u \). Then we can solve a linear PDE for the correction \( \delta u \) and use (23) to find a new approximation \( u^{n,k+1}=u^{n,k}+\delta u \) to \( u^{n} \).
Inserting (23) in (20) gives $$ \begin{equation} \frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^{n,k} + \delta u)\nabla (u^{n,k}+\delta u)) + f(u^{n,k}+\delta u)\tp \tag{24} \end{equation} $$ We can Taylor expand \( \dfc(u^{n,k} + \delta u) \) and \( f(u^{n,k}+\delta u) \): $$ \begin{align*} \dfc(u^{n,k} + \delta u) & = \dfc(u^{n,k}) + \frac{d\dfc}{du}(u^{n,k}) \delta u + \Oof{\delta u^2}\approx \dfc(u^{n,k}) + \dfc^{\prime}(u^{n,k})\delta u,\\ f(u^{n,k}+\delta u) &= f(u^{n,k}) + \frac{df}{du}(u^{n,k})\delta u + \Oof{\delta u^2}\approx f(u^{n,k}) + f^{\prime}(u^{n,k})\delta u\tp \end{align*} $$ Inserting the linear approximations of \( \dfc \) and \( f \) in (24) results in $$ \begin{align} \frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} &= \nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k}) + f(u^{n,k}) + \nonumber\\ &\qquad \nabla\cdot (\dfc(u^{n,k})\nabla \delta u) + \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k}) + \nonumber\\ &\qquad \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla \delta u) + f^{\prime}(u^{n,k})\delta u\tp \tag{25} \end{align} $$ The term \( \dfc^{\prime}(u^{n,k})\delta u\nabla \delta u \) is \( \Oof{\delta u^2} \) and therefore omitted. Reorganizing the equation gives a PDE for \( \delta u \) that we can write in short form as $$ \delta F(\delta u; u^{n,k}) = -F(u^{n,k}),$$ where $$ \begin{align} F(u^{n,k}) &= \frac{u^{n,k} - u^{n-1}}{\Delta t} - \nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k}) + f(u^{n,k}), \tag{26}\\ \delta F(\delta u; u^{n,k}) &= - \frac{1}{\Delta t}\delta u + \nabla\cdot (\dfc(u^{n,k})\nabla \delta u) + \nonumber\\ &\quad \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k}) + f^{\prime}(u^{n,k})\delta u\tp \end{align} $$ Note that \( \delta F \) is a linear function of \( \delta u \), and \( F \) contains only terms that are known, such that the PDE for \( \delta u \) is indeed linear.
We can rewrite the PDE for \( \delta u \) in a slightly different way too if we define \( u^{n,k} + \delta u \) as \( u^{n,k+1} \). $$ \begin{align} & \frac{u^{n,k+1} - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k+1}) + f(u^{n,k})\nonumber\\ &\qquad + \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k}) + f^{\prime}(u^{n,k})\delta u\tp \end{align} $$ Note that the first line is the same PDE as arise in the Picard iteration, while the remaining terms arise from the differentiations that are an inherent ingredient in Newton's method.
For coding we want to introduce \( u \) for \( u^n \), \( u^{-} \) for \( u^{n,k} \) and \( u^{(1)} \) for \( u^{n-1} \). The formulas for \( F \) and \( \delta F \) are then more clearly written as $$ \begin{align} F(u^{-}) &= \frac{u^{-} - u^{(1)}}{\Delta t} - \nabla\cdot (\dfc(u^{-})\nabla u^{-}) + f(u^{-}), \tag{27}\\ \delta F(\delta u; u^{-}) &= - \frac{1}{\Delta t}\delta u + \nabla\cdot (\dfc(u^{-})\nabla \delta u) + \nonumber\\ &\quad \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla u^{-}) + f^{\prime}(u^{-})\delta u\tp \end{align} $$ The form that orders the PDE as the Picard iteration terms plus the Newton method's derivative terms becomes $$ \begin{align} & \frac{u - u^{(1)}}{\Delta t} = \nabla\cdot (\dfc(u^{-})\nabla u) + f(u^{-}) + \nonumber\\ &\qquad \gamma(\nabla\cdot (\dfc^{\prime}(u^{-})(u - u^{-})\nabla u^{-}) + f^{\prime}(u^{-})(u - u^{-}))\tp \end{align} $$ The Picard and full Newton versions correspond to \( \gamma=0 \) and \( \gamma=1 \), respectively.
Some may prefer to derive the linearized PDE for \( \delta u \) using the more compact notation. We start with inserting \( u^n=u^{-}+\delta u \) to get $$ \frac{u^{-} +\delta u - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^{-} + \delta u)\nabla (u^{-}+\delta u)) + f(u^{-}+\delta u)\tp $$ Taylor expanding, $$ \begin{align*} \dfc(u^{-} + \delta u) & \approx \dfc(u^{-}) + \dfc^{\prime}(u^{-})\delta u,\\ f(u^{-}+\delta u) & \approx f(u^{-}) + f^{\prime}(u^{-})\delta u, \end{align*} $$ and inserting these expressions gives a less cluttered PDE for \( \delta u \): $$ \begin{align*} \frac{u^{-} +\delta u - u^{n-1}}{\Delta t} &= \nabla\cdot (\dfc(u^{-})\nabla u^{-}) + f(u^{-}) + \\ &\qquad \nabla\cdot (\dfc(u^{-})\nabla \delta u) + \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla u^{-}) + \\ &\qquad \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla \delta u) + f^{\prime}(u^{-})\delta u\tp \end{align*} $$
A Crank-Nicolson discretization of (17) applies a centered difference at \( t_{n+\frac{1}{2}} \): $$ [D_t u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^{n+\frac{1}{2}}\tp$$ Since \( u \) is not known at \( t_{n+\frac{1}{2}} \) we need to express the terms on the right-hand side via unknowns \( u^n \) and \( u^{n+1} \). The standard technique is to apply an arithmetic average, $$ u^{n+\frac{1}{2}}\approx \frac{1}{2}(u^n + u^{n+1})\tp$$ However, with nonlinear terms we have many choices of formulating an arithmetic mean: $$ \begin{align} [f(u)]^{n+\frac{1}{2}} &\approx f(\frac{1}{2}(u^n + u^{n+1})) = [f(\overline{u}^t)]^{n+\frac{1}{2}},\\ [f(u)]^{n+\frac{1}{2}} &\approx \frac{1}{2}(f(u^n) + f(u^{n+1})) =[\overline{f(u)}^t]^{n+\frac{1}{2}},\\ [\dfc(u)\nabla u]^{n+\frac{1}{2}} &\approx \dfc(\frac{1}{2}(u^n + u^{n+1}))\nabla (\frac{1}{2}(u^n + u^{n+1})) = [\dfc(\overline{u}^t)\nabla \overline{u}^t]^{n+\frac{1}{2}},\\ [\dfc(u)\nabla u]^{n+\frac{1}{2}} &\approx \frac{1}{2}(\dfc(u^n) + \dfc(u^{n+1}))\nabla (\frac{1}{2}(u^n + u^{n+1})) = [\overline{\dfc(u)}^t\nabla\overline{u}^t]^{n+\frac{1}{2}},\\ [\dfc(u)\nabla u]^{n+\frac{1}{2}} &\approx \frac{1}{2}(\dfc(u^n)\nabla u^n + \dfc(u^{n+1})\nabla u^{n+1}) = [\overline{\dfc(u)\nabla u}^t]^{n+\frac{1}{2}}\tp \end{align} $$
A big question is whether there are significant differences in accuracy between taking the products of arithmetic means or taking the arithmetic mean of products. Exercise 5: Find the truncation error of arithmetic mean of products investigates this question, and the answer is that the approximation is \( \Oof{\Delta t^2} \) in both cases.