This post builds on the formulations I showed in my previous posts by introducing their nonlinear versions.

In a typical nonlinear problem, the variational setting leads to the weak formulation

Find $u\in V$ such that

$$\begin{equation} F(u,v) = 0 \htmlId{eq:femnonlinear2}{} \tag{1}\end{equation}$$

for all $v\in V$ where the semilinear form $F$ is nonlinear in terms of $u$ and linear in terms of $v$.

We linearize $F$:

$$\begin{equation} \Lin [F(u,v)]_{u=\bar{u}} = F(\bar{u}, v) + \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}} \htmlId{eq:femnonlinear4}{} \tag{2}\end{equation}$$

Equating (2) to zero yields a linear system in terms of $\Var u$

$$\begin{equation} \boxed{ a(\Var u, v) = b(v) } \htmlId{eq:femnonlinear6}{} \tag{3}\end{equation}$$

where

$$\begin{equation} \boxed{ \begin{aligned} a(\Var u, v) &= \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}} \\ b(v) &= -F(\bar{u}, v) . \end{aligned} } \htmlId{eq:femnonlinear5}{} \tag{4}\end{equation}$$

We can compute the components of the matrices and vectors according to (3)

$$\begin{equation} \boxed{ \begin{alignedat}{3} \Aelid{I\!J}{} &= a(N^J,N^I) &&= \varn{F(u,N^I)}{u}{N^J}\evat_{u=\bar{u}} \\ b^{I} &= b(N^I) &&= -F(\bar{u}, N^I). \end{alignedat} } \htmlId{eq:femnonlinear9}{} \tag{5}\end{equation}$$

Then the update vector $\Var \Bu = [\Var u^1, \Var u^2, \dots, \Var u^{\nnode}]\tra$ is obtained by solving

$$\begin{equation} \BA \Var \Bu = \Bb \tag{6}\end{equation}$$

Letting $\Var u$ be the difference between consequent iterates, we obtain the update equation as

$$\begin{equation} \boxed{ \Bu \leftarrow \bar{\Bu} + \Var\Bu } \tag{7}\end{equation}$$

Example: Nonlinear Poisson’s Equation

Consider the following nonlinear Poisson’s equation

$$\begin{equation} \begin{alignedat}{4} - \nabla \dtp (g(u)\nabla u) &= f \quad && \text{in} \quad && \Omega \\ u &= 0 \quad && \text{on} \quad && \del\Omega \end{alignedat} \htmlId{eq:femnonlinear8}{} \tag{8}\end{equation}$$

The weak formulation reads

Find $u\in V$ such that

$$\begin{equation} - \int_\Omega \nabla \dtp (g(u)\nabla u) v \dv= \int_\Omega f v \dv \tag{9}\end{equation}$$

for all $v\in V$ where $V=H^1_0(\Omega)$.

Applying integration by parts and divergence theorem on the left-hand side

$$\begin{equation} \begin{aligned} \int_\Omega \nabla \dtp (g(u)\nabla u) v \dv &= \int_\Omega \nabla \dtp (g(u)\nabla (u) v) \dv - \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\ &= \underbrace{\int_{\del\Omega} g(u) v (\nabla u\dtp\Bn) \da}_{v = 0 \text{ on } \del\Omega} - \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\ \end{aligned} \tag{10}\end{equation}$$

Thus we have the semilinear form

$$\begin{equation} F(u,v) = \int_{\Omega} g(u) \nabla u \dtp \nabla v \dv - \int_{\Omega} f \, v \dv = 0 \tag{11}\end{equation}$$

The linearized version of this problem is then with (4)

$$\begin{equation} \begin{aligned} a(\Var u,v) &= \int_{\Omega} \rbr{\deriv{g}{u}\evat_{\bar{u}} \Var u\, \nabla \bar{u} + g(\bar{u})\nabla(\Var u)} \dtp \nabla v \dv \\ b(v) &= \int_{\Omega} [f \, v - g(\bar{u}) \nabla \bar{u} \dtp \nabla v] \dv\\ \end{aligned} \tag{12}\end{equation}$$

and the matrix and vector components are with (5)

$$\begin{equation} \begin{aligned} \Aelid{I\!J}{} &= \int_{\Omega} \rbr{\deriv{g}{u}\evat_{\bar{u}} N^J \, \nabla \bar{u} + g(\bar{u})\BB^J} \dtp \BB^I \dv \\ b^{I} &= \int_{\Omega} [f \, N^I - g(\bar{u}) \nabla \bar{u} \dtp \BB^I ] \dv\\ \end{aligned} \tag{13}\end{equation}$$

where the previous solution and its gradient are computed as

$$\begin{equation} \bar{u} = \suml{I=1}{\nnode} \bar{u}^I N^I \eqand \nabla \bar{u} = \suml{I=1}{\nnode} \bar{u}^I \BB^I . \tag{14}\end{equation}$$

Nonlinear Time-Dependent Problems

In the case of a nonlinear time-dependent problem, we have the following weak form:

Find $u \in V$ such that

$$\begin{equation} m(\dot{u}, v; t) + F(u,v; t) = 0 \htmlId{eq:nonlineartimedependentweak1}{} \tag{15}\end{equation}$$

for all $v \in V$ and $t \in [0,\infty)$ where $F$ is a semilinear form.

Discretization yields the following nonlinear system of equations

$$\begin{equation} \BM(t)\Bu + \Bf(u; t) = \Bzero \tag{16}\end{equation}$$

where

$$\begin{equation} \begin{aligned} M^{I\!J}(t) &= m(N^J, N^I; t) \\ f^{I}(u;t) &= F(u, N^I; t). \end{aligned} \tag{17}\end{equation}$$

Explicit Euler Scheme

We discretize in time with the finite difference $\dot{u} \approx [u_{n+1}-u_n]/{\Delta t}$ and linearity allows us to write

\begin{equation} \boxed{ m(\dot{u}, v; t) \approx \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] } \label{eq:discretetimedependent1} \end{equation}

We discretize the variational forms in time according to \eqref{eq:discretetimedependent1}, and evaluate the remaining terms at $t_n$:

$$\begin{equation} \frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_n, v; t_n) = 0 \tag{18}\end{equation}$$

The corresponding system of equations is

$$\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \Bf_n = \Bzero \tag{19}\end{equation}$$

where $\Bf_n = \Bf(u_n, t_n)$. This yields the following update equation

$$\begin{equation} \boxed{ \Bu_{n+1} = \BM_{n+1}\inv [\BM_n\Bu_n - \Delta t \Bf_n] } \tag{20}\end{equation}$$

For a time-independent $m$, this becomes

$$\begin{equation} \Bu_{n+1} = \Bu_n - \Delta t \BM\inv\Bf_n \tag{21}\end{equation}$$

Implicit Euler Scheme

For the implicit scheme, we evaluate the remaining terms at $t_{n+1}$ and let the result be equal to

$$\begin{equation} G(u_{n+1}, v) := \frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_{n+1}, v; t_{n+1}) = 0 \tag{22}\end{equation}$$

We will hereon replace $u_{n+1}$ with $u$ for brevity. The update of this nonlinear system requires the linearization of $G(u, v)$:

$$\begin{equation} \Lin[G(u,v)]_{u=\bar{u}} = G(\bar{u}, v) + \varn{G}{u}{\Var u}\evat_{u=\bar{u}} = 0 \tag{23}\end{equation}$$

We thus have the following linear setting for the Newton update $\Var u$:

$$\begin{equation} a(\Var u, v) = b(v) \tag{24}\end{equation}$$

where

$$\begin{equation} \begin{aligned} a(\Var u, v) &:= \varn{G}{u}{\Var u} \evat_{u=\bar{u}} = \frac{1}{\Delta t} m(\Var u, v; t_{n+1}) + \varn{F(u, v; t_{n+1})}{u}{\Var u} \evat_{u=\bar{u}} \\ b(v) &:= -G(\bar{u}, v) = - F(\bar{u}, v; t_{n+1}) -\frac{1}{\Delta t} [m(\bar{u},v;t_{n+1}) - m(u_{n},v;t_{n})] \end{aligned} \tag{25}\end{equation}$$

Discretization yields

$$\begin{equation} \rbr{\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}}\Var \Bu = \Bb \tag{26}\end{equation}$$

where

$$\begin{equation} \tilde{A}^{I\!J} := \varn{F(u, N^I;t_{n+1})}{u}{N^J} \evat_{u=\bar{u}} \eqand b^I := b(N^I) \tag{27}\end{equation}$$

The Newton update is rendered

$$\begin{equation} \boxed{ \Bu \leftarrow \bar{\Bu} + \Var\Bu \eqwith \Var \Bu = [\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}]\inv\Bb } \tag{28}\end{equation}$$

which is repeated until the solution for the next timestep $\Bu$ converges to a satisfactory value.

Nonlinear Coupled Problems

For a nonlinear coupled problem, the weak formulation is as follows

Find $u\in V_1$, $y\in V_2$ such that

$$\begin{equation} \begin{aligned} F(u, y, v) &= 0 \\ G(u, y, w) &= 0 \\ \end{aligned} \htmlId{eq:nonlinearcoupled1}{} \tag{29}\end{equation}$$

for all $v\in V_1$, $w \in V_2$ where $F(\cdot,\cdot, \cdot)$, $G(\cdot, \cdot, \cdot)$ are nonlinear in terms of $u$ and $y$ and linear in terms of $v$ and $w$.

We linearize the semilinear forms about the nonlinear terms:

$$\begin{equation} \begin{alignedat}{4} \Lin[F(u, y, v)]_{\bar{u},\bar{y}} &= F(\bar{u},\bar{y},v) &&+ \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}} &&+ \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \\ \Lin[G(u, y, w)]_{\bar{u},\bar{y}} &= G(\bar{u},\bar{y},w) &&+ \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}} &&+ \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \end{alignedat} \htmlId{eq:nonlinearcoupled2}{} \tag{30}\end{equation}$$

where the evaluations take place at $u=\bar{u}$ and $y=\bar{y}$.

Equating the linearized residuals to zero, we obtain a linear system of the form

$$\begin{equation} \begin{alignedat}{3} a(u, v) &+ b(y, v) &&= c(v) \\ d(u, w) &+ e(y, w) &&= f(w) \\ \end{alignedat} \htmlId{eq:coupledweakform1}{} \tag{31}\end{equation}$$

with the bilinear forms $a$, $b$, $d$, $e$ and the linear forms $c$, $f$ which are defined as

$$\begin{equation} \begin{gathered} \begin{alignedat}{4} a(\Var u, v) &:= \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}} \quad & b(\Var y, v) &:= \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \\ d(\Var u, w) &:= \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}} \quad & e(\Var y, w) &:= \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \end{alignedat} \\ \text{and} \\ \begin{aligned} c(v) &:= -F(\bar{u},\bar{y},v) \\ f(w) &:= -G(\bar{u}, \bar{y}, w) \end{aligned} \end{gathered} \tag{32}\end{equation}$$

Discretizing as done in the previous section, we obtain the following linear system of equations

$$\begin{equation} \begin{bmatrix} \BA & \BB \\ \BD & \BE \end{bmatrix} \begin{bmatrix} \Var \Bu \\ \Var \By \end{bmatrix} = \begin{bmatrix} \Bc \\ \Bf \end{bmatrix} \tag{33}\end{equation}$$

whose solution yields the update values $\Var \Bu$ and $\Var \By$. Thus the Newton update equations are

$$\begin{equation} \begin{alignedat}{3} \Bu &\leftarrow \bar{\Bu} &&+ \Var\Bu \\ \By &\leftarrow \bar{\By} &&+ \Var\By . \end{alignedat} \tag{34}\end{equation}$$

Example: Cahn-Hilliard Equation

The Cahn-Hilliard equation describes the process of phase separation, by which the two components of a binary fluid spontaneously separate and form domains pure in each component. The problem is nonlinear, coupled and time-dependent. The IBVP reads

$$\begin{equation} \begin{alignedat}{4} \partd{c}{t} &= \nabla\dtp(\BM\nabla \mu) \qquad&& \text{in} \qquad&& \Omega\times I \\ \nabla c\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\ \nabla \mu\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\ c &= c_0 && \text{in} && \Omega, t = 0 \\ \mu &= 0 && \text{in} && \Omega, t = 0 \\ \end{alignedat} \htmlId{eq:cahnhilliard1}{} \tag{35}\end{equation}$$

where

$$\begin{equation} \mu = \deriv{f}{c} - \nabla\dtp(\BLambda\nabla c) \htmlId{eq:cahnhilliard2}{} \tag{36}\end{equation}$$

and $t\in I = [0,\infty)$. Here,

• $c$ is the scalar variable for concentration,
• $\mu$ is the scalar variable for the chemical potential,
• $f: c \mapsto f(c)$ is the function representing chemical free energy,
• $\BM$ is a second-order tensor describing the mobility of the chemical,
• $\BLambda$ is a second-order tensor describing both the interface thickness and direction of phase transition.

The fourth-order PDE governing the problem can be formulated as a coupled system of two second-order PDEs with the variables $c$ and $\mu$, as demonstrated in (35) and (36).

The weak formulation then reads

Find $c \in V_1$, $\mu\in V_2$ such that

$$\begin{equation} \begin{aligned} \int_\Omega \partd{c}{t} v \dx - \int_\Omega \nabla\dtp(\BM\nabla \mu) v \dx &=0 \\ \int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx + \int_\Omega \nabla\dtp(\BLambda\nabla c) w\dx &= 0 \end{aligned} \tag{37}\end{equation}$$

for all $v \in V_1$, $w \in V_2$ and $t \in I$.

We discretize in time implicitly with $\del c/\del t \approx (c_{n+1}-c_n)/\Var t$. We also denote the values for the next timestep $c_{n+1}$ and $\mu_{n+1}$ as $c$ and $\mu$ for brevity. Using integration-by-parts, the divergence theorem, and the given boundary conditions, we arrive at the following nonlinear forms

$$\begin{equation} \begin{alignedat}{3} F(c,\mu,v) &= \int_\Omega \frac{1}{\Var t} (c-c_n) v \dx + \int_\Omega (\BM\nabla \mu)\dtp \nabla v \dx &&= 0 \\ G(c,\mu,w) &= \int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx - \int_\Omega (\BLambda\nabla c)\dtp \nabla w\dx &&= 0 \end{alignedat} \tag{38}\end{equation}$$

which is a nonlinear coupled system of the form (29).

We linearize the forms according to (30) and obtain the following variations

$$\begin{align*} \varn{F}{c}{\Var c} &= \int_\Omega \frac{1}{\Var t} \Var c\, v \dx \\ \varn{F}{\mu}{\Var \mu} &= \int_\Omega (\BM\nabla (\Var\mu))\dtp \nabla v \dx \\ \varn{G}{c}{\Var c} &= - \int_\Omega \dderiv{f}{c}\Var c \, w \dx - \int_\Omega (\BLambda\nabla (\Var c))\dtp \nabla w\dx \\ \varn{G}{\mu}{\Var \mu} &= \int_\Omega \Var\mu \, w \dx \end{align*}$$

We substitute basis functions and obtain our system matrix and vectors

$$\begin{align*} P^{I\!J} &= \int_\Omega \frac{1}{\Var t} N^JN^I \dx \\ Q^{IL} &= \int_\Omega (\BM\BB^L)\dtp\BB^I \dx \\ r^{I} &= \int_\Omega \frac{1}{\Var t}(\bar{c}-c_n)N^I \dx + \int_\Omega (\BM\nabla\bar{\mu})\dtp\BB^I\dx \\ S^{K\!J} &= - \int_\Omega \dderiv{f}{c}\evat_{c=\bar{c}} N^J N^K \dx - \int_\Omega (\BLambda \BB^J)\dtp \BB^K\dx \\ T^{K\!L} &= \int_\Omega N^L N^K \dx \\ u^{K} &= \int_\Omega \sbr{\bar{\mu} - \deriv{f}{c}\evat_{c=\bar{c}}} N^K \dx - \int_\Omega (\BLambda\nabla \bar{c})\dtp \BB^K\dx \end{align*}$$

which constitute the system

$$\begin{equation} \begin{bmatrix} \BP & \BQ \\ \BS & \BT \end{bmatrix} \begin{bmatrix} \Var \Bc \\ \Var \Bmu \end{bmatrix} = \begin{bmatrix} \Br \\ \Bu \end{bmatrix} \tag{39}\end{equation}$$

Solution yields the update values $\Var \Bc$ and $\Var \Bmu$. The Newton update equations are then

$$\begin{equation} \begin{alignedat}{3} \Bc &\leftarrow \bar{\Bc} &&+ \Var\Bc \\ \Bmu &\leftarrow \bar{\Bmu} &&+ \Var\Bmu . \end{alignedat} \tag{40}\end{equation}$$

The system is solved for $c_{n+1}$ and $\mu_{n+1}$ at each $t=t_n$ to obtain the evolutions of the concentration and chemical potential.

In this post, I’ll introduce the FE formulation of a generalized linear and coupled weak form. Said weak formulation has the form

Find $u\in V_1$, $y\in V_2$ such that

$$\begin{equation} \begin{alignedat}{3} a(u, v) &+ b(y, v) &&= c(v) \\ d(u, w) &+ e(y, w) &&= f(w) \\ \end{alignedat} \htmlId{eq:coupledweakform1}{} \tag{1}\end{equation}$$

for all $v\in V_1$, $w \in V_2$ where $a(\cdot, \cdot): V_1\times V_1 \to \IR$, $b(\cdot, \cdot): V_2\times V_1 \to \IR$, $d(\cdot, \cdot): V_1\times V_2 \to \IR$, $e(\cdot, \cdot): V_2\times V_2 \to \IR$ are bilinear forms and $c(\cdot): V_1\to \IR$, $f(\cdot): V_2\to \IR$ are linear forms.

Here, the objective is to solve for the two unknown functions $u$ and $y$. One can also imagine an arbitrary degree of coupling between $n$ variables with $n$ equations.

We introduce the following discretizations

$$\begin{equation} \begin{alignedat}{3} u_h &= \suml{J=1}{n_n^1} u^J N^J \qquad\qquad & v_h &= \suml{I=1}{n_n^1} u^I N^I \qquad\qquad & u_h, v_h\in V_{h1} \\ y_h &= \suml{L=1}{n_n^2} u^L N^L & w_h &= \suml{K=1}{n_n^2} u^K N^K & y_h, w_h\in V_{h2} \\ \end{alignedat} \tag{2}\end{equation}$$

where the corresponding number of shape functions are $n_n^1$ and $n_n^2$, respectively.

Substituting the discretizations in (1), we obtain two linear systems of equations

$$\begin{equation} \begin{alignedat}{3} \suml{J=1}{n_n^1} a(N^J, N^I) \,u^J &+ \suml{L=1}{n_n^2} b(N^L, N^I) \, y^L &&= c(N^I) \\ \suml{J=1}{n_n^1} d(N^J, N^K) \,u^J &+ \suml{L=1}{n_n^2} e(N^L, N^K) \, y^L &&= f(N^K) \\ \end{alignedat} \tag{3}\end{equation}$$

for $I=1,\dots,n_n^1$ and $K=1,\dots,n_n^2$.

We write this system as

$$\begin{equation} \boxed{ \begin{alignedat}{3} \BA \Bu &+ \BB\By &&= \Bc \\ \BD \Bu &+ \BE\By &&= \Bf \\ \end{alignedat} \eqor \begin{bmatrix} \BA & \BB \\ \BD & \BE \end{bmatrix} \begin{bmatrix} \Bu \\ \By \end{bmatrix} = \begin{bmatrix} \Bc \\ \Bf \end{bmatrix} } \htmlId{eq:coupledsystem1}{} \tag{4}\end{equation}$$

where the components of given matrices and vectors are defined as

$$\begin{equation} \begin{alignedat}{6} A^{I\!J} &:= a(N^J, N^I) \qquad & B^{I\!L} &:= b(N^L, N^I) \qquad & c^{I} &:= c(N^I)\\ D^{K\!J} &:= d(N^J, N^K) & E^{K\!L} &:= e(N^L, N^K) & f^{K} &:= f(N^K)\\ \end{alignedat} \tag{5}\end{equation}$$

Solution of (4) yields the unknown vectors $\Bu$ and $\By$.

Time dependent problems are commonplace in physics, chemistry and many other disciplines. In this post, I’ll introduce the FE formulation of linear time-dependent problems and derive formulas for explicit and implicit Euler integration.

The weak formulation of a first order time-dependent problem reads:

Find $u \in V$ such that

$$\begin{equation} m(\dot{u}, v; t) + a(u,v; t) = b(v; t) \htmlId{eq:timedependentweak1}{} \tag{1}\end{equation}$$

for all $v \in V$ and $t \in [0,\infty)$.

We can convert (1) into a system of equations

$$\begin{equation} \BM(t)\dot{\Bu} + \BA(t)\Bu = \Bb(t) \tag{2}\end{equation}$$

where the components of the matrices and vectors involved are calculated as

$$\begin{equation} \begin{aligned} M^{I\!J}(t) &= m(N^J, N^I; t) \\ A^{I\!J}(t) &= a(N^J, N^I; t) \\ b^{I}(t) &= b(N^I; t). \end{aligned} \tag{3}\end{equation}$$

If we further discretize in time with the finite difference $\dot{u} \approx [u_{n+1}-u_n]/{\Delta t}$, linearity allows us to write

$$\begin{equation} \boxed{ m(\dot{u}, v; t) \approx \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] } \htmlId{eq:discretetimedependent1}{} \tag{4}\end{equation}$$

This reflects on the system as

$$\begin{equation} \BM(t)\dot{\Bu} \approx \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] \htmlId{eq:discretetimedependent2}{} \tag{5}\end{equation}$$

Here, $u_{n+1}:= u(x, t_{n+1})$, $\BM_{n+1} = \BM(t_{n+1})$ and vice versa for $u_n$ and $\BM_n$.

Explicit Euler Scheme

For the explicit Euler scheme, we substitute evaluate the remaining terms at $t_n$

$$\begin{equation} \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] + a(u_n,v; t_n) = b(v; t_n) \quad \forall v \in V\,. \tag{6}\end{equation}$$

The corresponding system is

$$\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \BA_n\Bu_n = \Bb_n \tag{7}\end{equation}$$

The update equation becomes

$$\begin{equation} \boxed{ \Bu_{n+1} = \BM_{n+1}\inv [\BM_n\Bu_n + \Delta t(\Bb_n - \BA_n\Bu_n)] } \tag{8}\end{equation}$$

If $m$ is time-independent, that is $m(\dot{u}, v;t) = m(\dot{u}, v)$, we have

$$\begin{equation} \Bu_{n+1} = \Bu_n + \Delta t\, \BM\inv(\Bb_n - \BA_n\Bu_n) \tag{9}\end{equation}$$

Implicit Euler Scheme

For the implicit Euler scheme, we substitute evaluate the remaining terms at $t_{n+1}$

$$\begin{equation} \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] + a(u_n,v; t_{n+1}) = b(v; t_{n+1}) \quad \forall v \in V\,. \tag{10}\end{equation}$$

The corresponding system is

$$\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \BA_{n+1}\Bu_{n+1} = \Bb_{n+1} \tag{11}\end{equation}$$

The update equation becomes

$$\begin{equation} \boxed{ \Bu_{n+1} = [\BM_{n+1}+\Delta t \BA_{n+1}]\inv [\BM_n\Bu_n + \Delta t \,\Bb_{n+1}] } \tag{12}\end{equation}$$

If $m$ is time-independent, one can just substitute $\BM=\BM_{n+1}=\BM_n$.

Example: Reaction-Advection-Diffusion Equation

The IBVP of a linear reaction-advection-diffusion problem reads

$$\begin{equation} \begin{alignedat}{4} \partd{u}{t} &= \nabla\dtp(\BD\nabla u) - \nabla\dtp(\Bc u) + ru + f \qquad&& \text{in} \qquad&& \Omega\times I\\ u &= \bar{u} && \text{on} && \del\Omega\times I\\ u &= u_0 && \text{in} && \Omega, t = 0 \\ \end{alignedat} \tag{13}\end{equation}$$

where $t\in I = [0,\infty)$,

• $\BD$ is a second-order tensor describing the diffusivity of $u$,
• $\Bc$ is a vector describing the velocity of advection,
• $r$ is a scalar describing the rate of reaction,
• and $f$ is a source term for $u$.

The weak formulation is then

Find $u \in V$ such that

$$\begin{equation} \int_\Omega \dot{u} v \dv = \int_\Omega [\nabla\dtp(\BD\nabla u) - \nabla\dtp(\Bc u) + ru + f] v \dv \tag{14}\end{equation}$$

for all $v \in V$ and $t \in I$.

We have the following integration by parts relationships:

$$\begin{equation} \int_\Omega \nabla \dtp(\BD\nabla u) v \dv = \cancel{\int_\Omega \nabla\dtp(v\BD\nabla u) \dv} - \int_\Omega (\BD\nabla u)\dtp\nabla v \dv \end{equation}$$

for the diffusive part and

$$\begin{equation} \int_\Omega \nabla\dtp(\Bc u) v \dv = \cancel{\int_\Omega \nabla \dtp (\Bc u v) \dv} - \int_\Omega u \Bc \dtp \nabla v \dv \tag{15}\end{equation}$$

for the advective part. The canceled terms are due to divergence theorem and the fact that $v=0$ on the boundary. Then our variational formulation is of the form (1) where

$$\begin{align*} m(\dot{u}, v) &= \int_\Omega \dot{u} v \dv \\ a(u, v) &= \int_\Omega (\BD\nabla u) \dtp \nabla v \dv - \int_\Omega u\Bc \dtp \nabla v \dv - \int_\Omega ruv \dv \\ b(v) &= \int_\Omega fv \dv \end{align*}$$

From these forms, we obtain the following system matrices and vector

$$\begin{align*} M^{I\!J} &= \int_\Omega N^J N^I \dv \\ A^{I\!J} &= \int_\Omega (\BD\BB^J) \dtp \BB^I \dv - \int_\Omega N^J\Bc \dtp \BB^I \dv - \int_\Omega r N^JN^I \dv \\ b^I &= \int_\Omega f N^I \dv \end{align*}$$

where $\BM$ is constant through time.