Linearization at the differential equation level

The attention is now turned to nonlinear partial differential equations (PDEs) and application of the techniques explained for ODEs. The model problem is a nonlinear diffusion equation

(1)\[ \frac{\partial u}{\partial t} = \nabla\cdot ({\alpha}(u)\nabla u) + f(u),\quad \boldsymbol{x}\in\Omega,\ t\in (0,T],\]

(2)\[ -{\alpha}(u)\frac{\partial u}{\partial n} = g,\quad \boldsymbol{x}\in\partial\Omega_N,\ t\in (0,T],\]

(3)\[ u = u_0,\quad \boldsymbol{x}\in\partial\Omega_D,\ t\in (0,T]{\thinspace .}\]

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 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.

Explicit time integration

The nonlinearities in the PDE are trivial to deal with if we choose an explicit time integration method for (1), such as the Forward Euler method:

\[D_t^+ u = \nabla\cdot ({\alpha}(u)\nabla u) + f(u)]^n,\]

which leads to a linear equation in the unknown \(u^{n+1}\):

\[\frac{u^{n+1} - u^n}{\Delta t} = \nabla\cdot ({\alpha}(u^n)\nabla u^n) + f(u^n){\thinspace .}\]

Picard iteration (3)

A Backward Euler scheme for (1) reads

\[D_t^- u = \nabla\cdot ({\alpha}(u)\nabla u) + f(u)]^n{\thinspace .}\]

Written out,

(4)\[ \frac{u^{n} - u^{n-1}}{\Delta t} = \nabla\cdot ({\alpha}(u^n)\nabla u^n) + f(u^n)\]

This is a nonlinear, stationary PDE for the unknown function \(u^n(\boldsymbol{x})\). We introduce a Picard iteration with \(k\) as iteration counter. A typical linearization of the \(\nabla\cdot{\alpha}(u^n)\nabla u^n\) term in iteration \(k+1\) is to use the previously computed \(u^{n,k}\) approximation in the diffusion coefficient: \({\alpha}(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

(5)\[ \frac{u^{n,k+1} - u^{n-1}}{\Delta t} = \nabla\cdot ({\alpha}(u^{n,k}) \nabla u^{n,k+1}) + f(u^{n,k}){\thinspace .}\]

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 notation where \(u\) corresponds to the unknown we want to solve for, i.e., \(u^{n,k+1}\) above, let \(u_{-}\) be the most recently computed value, \(u^{n,k}\) above, and let \(u_1\) denote the unknown function at the previous time level, \(u^{n-1}\) above. The PDE to be solved in a Picard iteration then looks like

(6)\[ \frac{u - u_1}{\Delta t} = \nabla\cdot ({\alpha}(u_{-}) \nabla u) + f(u_{-}){\thinspace .}\]

At the beginning of the iteration we start with the value from the previous time level: \(u_{-}=u_1\).

Newton’s method (4)

At time level \(n\) we have to solve the stationary PDE (4), 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 \(u^n\). We seek a better approximation on the form

(7)\[ u^{n} = u^{n,k} + \delta u{\thinspace .}\]

The idea is to insert (7) in (4), Taylor expand the nonlinearities and only keep the terms that are linear in \(\delta u\). Then we can solve a linear PDE for the correction \(\delta u\) and use (7) to find a new approximation \(u^{n,k+1}=u^{n,k}+\delta u\) to \(u^{n}\).

Inserting (7) in (4) gives

(8)\[ \frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} = \nabla\cdot ({\alpha}(u^{n,k} + \delta u)\nabla (u^{n,k}+\delta u)) + f(u^{n,k}+\delta u)\]

We can Taylor expand \({\alpha}(u^{n,k} + \delta u)\) and \(f(u^{n,k}+\delta u)\):

\[\begin{split}{\alpha}(u^{n,k} + \delta u) & = {\alpha}(u^{n,k}) + \frac{d{\alpha}}{du}(u^{n,k}) \delta u + {\mathcal{O}(\delta u^2)}\approx {\alpha}(u^{n,k}) + {\alpha}'(u^{n,k})\delta u,\\ f(u^{n,k}+\delta u) &= f(u^{n,k}) + \frac{df}{du}(u^{n,k})\delta u + {\mathcal{O}(\delta u^2)}\approx f(u^{n,k}) + f'(u^{n,k})\delta u{\thinspace .}\end{split}\]

Inserting the linear approximations of \({\alpha}\) and \(f\) in (8) results in

\[\frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} = \nabla\cdot ({\alpha}(u^{n,k})\nabla u^{n,k}) + f(u^{m,k}) + \nonumber\]

\[\quad \nabla\cdot ({\alpha}(u^{n,k})\nabla \delta u) + \nabla\cdot ({\alpha}'(u^{n,k})\delta u\nabla u^{n,k}) + \nonumber\]

(9)\[ \quad \nabla\cdot ({\alpha}'(u^{n,k})\delta u\nabla \delta u) + f'(u^{n,k})\delta u\]

The term \({\alpha}'(u^{n,k})\delta u\nabla \delta u\) is \({\mathcal{O}(\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

(10)\[ F(u^{n,k}) = \frac{u^{n,k} - u^{n-1}}{\Delta t} - \nabla\cdot ({\alpha}(u^{n,k})\nabla u^{n,k}) + f(u^{n,k}),\]

\[\delta F(\delta u; u^{n,k}) = - \frac{1}{\Delta t}\delta u + \nabla\cdot ({\alpha}(u^{n,k})\nabla \delta u) + \nonumber\]

\[\quad \nabla\cdot ({\alpha}'(u^{n,k})\delta u\nabla u^{n,k}) + f'(u^{n,k})\delta u{\thinspace .}\]

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.

The form \(\delta F = -F\) resembles the Newton system \(J\delta u =-F\) for systems of algebraic equations, with \(\delta F\) as \(J\delta u\). The unknown vector in a linear system of algebraic equations enters the system as a matrix-vector product (\(J\delta u\)), while at the PDE level we have a linear differential operator instead (\(\delta F\)).

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}\).

\[ \frac{u^{n,k+1} - u^{n-1}}{\Delta t} = \nabla\cdot ({\alpha}(u^{n,k})\nabla u^{n,k+1}) + f(u^{n,k}) + \nonumber\]

\[\qquad \nabla\cdot ({\alpha}'(u^{n,k})\delta u\nabla u^{n,k}) + f'(u^{n,k})\delta u{\thinspace .}\]

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,k}\) and \(u_1\) for \(u^{n-1}\). The formulas for \(F\) and \(\delta F\) are then

(11)\[ F(u_{-}) = \frac{u_{-} - u_1}{\Delta t} - \nabla\cdot ({\alpha}(u_{-})\nabla u_{-}) + f(u_{-}),\]

\[\delta F(\delta u; u_{-}) = - \frac{1}{\Delta t}\delta u + \nabla\cdot ({\alpha}(u_{-})\nabla \delta u) + \nonumber\]

\[\quad \nabla\cdot ({\alpha}'(u_{-})\delta u\nabla u_{-}) + f'(u_{-})\delta u{\thinspace .}\]

The form that orders the PDE as the Picard iteration terms plus the Newton method’s derivative terms becomes

\[ \frac{u - u_1}{\Delta t} = \nabla\cdot ({\alpha}(u_{-})\nabla u) + f(u_{-}) + \nonumber\]

\[\qquad \nabla\cdot ({\alpha}'(u_{-})\delta u\nabla u_{-}) + f'(u_{-})\delta u{\thinspace .}\]