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.