Multi-dimensional PDE problems

Finite element discretization

The derivation of \(F_i\) and \(J_{i,j}\) in the 1D model problem is easily generalized to multi-dimensional problems. For example, Backward Euler discretization of the PDE

(1)\[ u_t = \nabla\cdot({\alpha}(u)\nabla u) + f(u),\]

gives the nonlinear time-discrete PDEs

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

or with \(u^n\) simply as \(u\) and \(u^{n-1}\) as \(u_1\),

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

The variational form, assuming homogeneous Neumann conditions for simplicity, becomes

(2)\[ \int_\Omega (uv + \Delta t{\alpha}(u)\nabla u\cdot\nabla v - \Delta t f(u)v - u_1v){\, \mathrm{d}x}{\thinspace .}\]

The nonlinear algebraic equations follow from setting \(v={\psi}_i\) and using the representation \(u=\sum_kc_k{\psi}_k\), which we just write as

(3)\[ F_i = \int_\Omega (u{\psi}_i + \Delta t{\alpha}(u)\nabla u\cdot\nabla {\psi}_i - \Delta t f(u){\psi}_i - u_1{\psi}_i){\, \mathrm{d}x}{\thinspace .}\]

Picard iteration needs a linearization where we use the most recent approximation \(u_{-}\) to \(u\) in \({\alpha}\) and \(f\):

(4)\[ F_i \approx \hat F_i = \int_\Omega (u{\psi}_i + \Delta t{\alpha}(u_{-})\nabla u\cdot\nabla {\psi}_i - \Delta t f(u_{-}){\psi}_i - u_1{\psi}_i){\, \mathrm{d}x}{\thinspace .}\]

The equations \(\hat F_i=0\) are now linear and we can easily derive a linear system for the unknown coefficients \(\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}\) by inserting \(u=\sum_jc_j{\psi}_j\).

In Newton’s method we need to evaluate \(F_i\) with the known value \(u_{-}\) for \(u\):

(5)\[ F_i \approx \hat F_i = \int_\Omega (u_{-}{\psi}_i + \Delta t{\alpha}(u_{-}) \nabla u_{-}\cdot\nabla {\psi}_i - \Delta t f(u_{-}){\psi}_i - u_1{\psi}_i){\, \mathrm{d}x}{\thinspace .}\]

The Jacobian is obtained by differentiating (3) and using \(\partial u/\partial c_j={\psi}_j\):

\[J_{i,j} = \frac{\partial F_i}{\partial c_j} = \int_\Omega ({\psi}_j{\psi}_i + \Delta t{\alpha}'(u){\psi}_j \nabla u\cdot\nabla {\psi}_i + \Delta t{\alpha}(u)\nabla{\psi}_j\cdot\nabla{\psi}_i - \nonumber\]

(6)\[ \ \Delta t f'(u){\psi}_j{\psi}_i){\, \mathrm{d}x}{\thinspace .}\]

The evaluation of \(J_{i,j}\) as the coefficient matrix in the linear system in Newton’s method applies the known approximation \(u_{-}\) for \(u\):

\[J_{i,j} = \int_\Omega ({\psi}_j{\psi}_i + \Delta t{\alpha}'(u_{-}){\psi}_j \nabla u_{-}\cdot\nabla {\psi}_i + \Delta t{\alpha}(u_{-})\nabla{\psi}_j\cdot\nabla{\psi}_i - \nonumber\]

(7)\[ \ \Delta t f'(u_{-}){\psi}_j{\psi}_i){\, \mathrm{d}x}{\thinspace .}\]

Hopefully, these example also show how convenient the notation with \(u\) and \(u_{-}\) is: the unknown to be computed is always \(u\) and linearization by inserting known (previously computed) values is a matter of adding an underscore. One can take great advantage of this quick notation in software [Ref2].

Non-homogeneous Neumann conditions

A natural physical flux condition for the PDE (1) takes the form of a non-homogeneous Neumann condition

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

where \(g\) is a prescribed function and \(\partial\Omega_N\) is a part of the boundary of the domain \(\Omega\). From integrating \(\int_\Omega\nabla\cdot({\alpha}\nabla u){\, \mathrm{d}x}\) by parts, we get a boundary term

(9)\[ \int_{\partial\Omega_N}{\alpha}(u)\frac{\partial u}{\partial u}v{\, \mathrm{d}s}{\thinspace .}\]

Inserting the condition (8) into this term results in an integral over prescribed values: \(-\int_{\partial\Omega_N}gv{\, \mathrm{d}s}\). The nonlinearity in the \({\alpha}(u)\) coefficient condition (8) therefore does not contribute with a nonlinearity in the variational form.

Robin conditions

Heat conduction problems often apply a kind of Newton’s cooling law, also known as a Robin condition, at the boundary:

(10)\[ -{\alpha}(u)\frac{\partial u}{\partial u} = h_T(u)(u-T_s(t)),\quad\boldsymbol{x}\in\partial\Omega_R,\]

where \(h_T(u)\) is a heat transfer coefficient between the body (\(\Omega\)) and its surroundings, \(T_s\) is the temperature of the surroundings, and \(\partial\Omega_R\) is a part of the boundary where this Robin condition applies. The boundary integral (9) now becomes

\[\int_{\partial\Omega_R}h_T(u)(u-T_s(T))v{\, \mathrm{d}s},\]

by replacing \({\alpha}(u)\partial u/\partial u\) by \(h_T(u-T_s)\). Often, \(h_T(u)\) can be taken as constant, and then the boundary term is linear in \(u\), otherwise it is nonlinear and contributes to the Jacobian in a Newton method. Linearization in a Picard method will typically use a known value in \(h_T\), but keep the \(u\) in \(u-T_s\) as unknown: \(h_T(u_{-})(u-T_s(t))\).

Finite difference discretization

A typical diffusion equation

\[u_t = \nabla\cdot({\alpha}(u)\nabla u) + f(u),\]

can be discretized by (e.g.) a Backward Euler scheme, which in 2D can be written

\[[D_t^- u = D_x\overline{{\alpha}}^xD_x u + D_y\overline{{\alpha}}^yD_y u + f(u)]_{i,j}^n{\thinspace .}\]

We do not dive into details of boundary conditions now. Dirichlet and Neumann conditions are handled as in linear diffusion problems.

Writing the scheme out, putting the unknown values on the left-hand side and known values on the right-hand side, and introducing \(\Delta x=\Delta y=h\) to save some writing, one gets

\[\begin{split}u^n_{i,j} &- \frac{\Delta t}{h^2}( \frac{1}{2}({\alpha}(u_{i,j}^n) + {\alpha}(u_{i+1,j}^n))(u_{i+1,j}^n-u_{i,j}^n)\\ &\quad - \frac{1}{2}({\alpha}(u_{i-1,j}^n) + {\alpha}(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j}^n) \\ &\quad + \frac{1}{2}({\alpha}(u_{i,j}^n) + {\alpha}(u_{i,j+1}^n))(u_{i,j+1}^n-u_{i,j}^n)\\ &\quad - \frac{1}{2}({\alpha}(u_{i,j-1}^n) + {\alpha}(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j-1}^n)) - \Delta tf(u_{i,j}^n) = u^{n-1}_{i,j}\end{split}\]

This defines a nonlinear algebraic system \(A(u)u=b(u)\). A Picard iteration applies old values \(u_{-}\) in \({\alpha}\) and \(f\), or equivalently, old values for \(u\) in \(A(u)\) and \(b(u)\). The result is a linear system of the same type as those arising from \(u_t = \nabla\cdot ({\alpha}(\boldsymbol{x})\nabla u) + f(\boldsymbol{x},t)\).

Newton’s method is as usual more involved. We first define the nonlinear algebraic equations to be solved, drop the superscript \(n\), and introduce \(u_1\) for \(u^{n-1}\):

\[\begin{split}F_{i,j} &= u^n_{i,j} - \frac{\Delta t}{h^2}(\\ &\quad \frac{1}{2}({\alpha}(u_{i,j}^n) + {\alpha}(u_{i+1,j}^n))(u_{i+1,j}^n-u_{i,j}^n) - \frac{1}{2}({\alpha}(u_{i-1,j}^n) + {\alpha}(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j}^n) + \\ &\quad \frac{1}{2}({\alpha}(u_{i,j}^n) + {\alpha}(u_{i,j+1}^n))(u_{i,j+1}^n-u_{i,j}^n) - \frac{1}{2}({\alpha}(u_{i,j-1}^n) + {\alpha}(u_{i,j}^n))(u_{i,j}^n-u_{i-1,j-1}^n)) -\\ &\quad \Delta tf(u_{i,j}^n) - u^{n-1}_{i,j} = 0{\thinspace .}\end{split}\]

It is convenient to work with two indices \(i\) and \(j\) in 2D finite difference discretizations, but it complicates the derivation of the Jacobian, which then gets four indices. The left-hand expression of an equation \(F_{i,j}=0\) is to be differentiated with respect to each of the unknowns \(u_{r,s}\) (short for \(u_{r,s}^n\)), \(r\in{\mathcal{I}_x}\), \(s\in{\mathcal{I}_y}\),

\[J_{i,j,r,s} = \frac{\partial F_{i,j}}{\partial u_{r,s}}{\thinspace .}\]

Given \(i\) and \(j\), only a few \(r\) and \(s\) indices give nonzero contribution since \(F_{i,j}\) contains \(u_{i\pm 1,j}\), \(u_{i,j\pm 1}\), and \(u_{i,j}\). Therefore, \(J_{i,j,r,s}\) is very sparse and we can set up the left-hand side of the Newton system as

\[\begin{split} J_{i,j,r,s}\delta u_{r,s} = J_{i,j,i,j}\delta u_{i,j} & + J_{i,j,i-1,j}\delta u_{i-1,j} + J_{i,j,i+1,j}\delta u_{i+1,j} + J_{i,j,i,j-1}\delta u_{i,j-1}\\ & + J_{i,j,i,j+1}\delta u_{i,j+1}\end{split}\]

The specific derivatives become

\[\begin{split}J_{i,j,i-1,j} &= \frac{\partial F_{i,j}}{\partial u_{i-1,j}}\\ &= \frac{\Delta t}{h^2}({\alpha}'(u_{i-1,j})(u_{i,j}-u_{i-1,j}) + {\alpha}(u_{i-1,j})(-1))\\ J_{i,j,i+1,j} &= \frac{\partial F_{i,j}}{\partial u_{i+1,j}}\\ &= \frac{\Delta t}{h^2}(-{\alpha}'(u_{i+1,j})(u_{i+1,j}-u_{i,j}) - {\alpha}(u_{i-1,j}))\\ J_{i,j,i,j-1} &= \frac{\partial F_{i,j}}{\partial u_{i,j-1}}\\ &= \frac{\Delta t}{h^2}({\alpha}'(u_{i,j-1})(u_{i,j}-u_{i,j-1}) + {\alpha}(u_{i,j-1})(-1))\\ J_{i,j,i,j+1} &= \frac{\partial F_{i,j}}{\partial u_{i,j+1}}\\ &= \frac{\Delta t}{h^2}(-{\alpha}'(u_{i,j+1})(u_{i,j+1}-u_{i,j}) - {\alpha}(u_{i,j-1}))\end{split}\]

The \(J_{i,j,i,j}\) entry has a few more terms. Inserting \(u_{-}\) for \(u\) in the \(J\) formula and then forming \(J\delta u=-F\) gives the linear system to be solved in each Newton iteration.

Continuation methods

Picard iteration or Newton’s method may diverge when solving PDEs with severe nonlinearities. Relaxation with \(\omega <1\) may help, but in highly nonlinear problems it can be necessary to introduce a continuation parameter \(\Lambda\) in the problem: \(\Lambda =0\) gives a version of the problem that is easy to solve, while \(\Lambda =1\) is the target problem. The idea is then to increase \(\Lambda\) in steps, \(\Lambda_0=0 ,\Lambda_1 <\cdots <\Lambda_n=1\), and use the solution from the problem with \(\Lambda_{i-1}\) as initial guess for the iterations in the problem corresponding to \(\Lambda_i\).

The continuation method is easiest to understand through an example. Suppose we intend to solve

\[-\nabla\cdot\left( ||\nabla u||^q\nabla u\right) = f,\]

which is an equation modeling the flow of a non-Newtonian fluid through i channel or pipe. For \(q=0\) we have the Poisson equation (corresponding to a Newtonian fluid) and the problem is linear. A typical value for pseudo-plastic fluids may be \(q_n=-0.8\). We can introduce the continuation parameter \(\Lambda\in [0,1]\) such that \(q=q_n\Lambda\). Let \(\{\Lambda_\ell\}_{\ell=0}^n\) be the sequence of \(\Lambda\) values in \([0,1]\), with corresponding \(q\) values \(\{q_\ell\}_{\ell=0}^n\). We can then solve a sequence of problems

\[-\nabla\cdot\left( ||\nabla u||^q_\ell\nabla u^\ell\right) = f,\quad \ell = 0,\ldots,n,\]

where the initial guess for iterating on \(u^{\ell}\) is the previously computed solution \(u^{\ell-1}\). If a particular \(\Lambda_\ell\) leads to convergence problems, one may try a smaller increase in \(\Lambda\): \(\Lambda_* = \frac{1}{2} (\Lambda_{\ell-1}+\Lambda_\ell)\), and repeat halving the step in \(\Lambda\) until convergence is reestablished.