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

$$$m(\dot{u}, v; t) + a(u,v; t) = b(v; t) \label{eq:timedependentweak1}$$$

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

We can convert \eqref{eq:timedependentweak1} into a system of equations

$$$\BM(t)\dot{\Bu} + \BA(t)\Bu = \Bb(t)$$$

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

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

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

$$$\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}$$$

This reflects on the system as

$$$\BM(t)\dot{\Bu} \approx \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] \label{eq:discretetimedependent2}$$$

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$

$$$\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\,.$$$

The corresponding system is

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

The update equation becomes

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

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

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

## Implicit Euler Scheme

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

$$$\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\,.$$$

The corresponding system is

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

The update equation becomes

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

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

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

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

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

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

We have the following integration by parts relationships:

$\require{cancel}$$\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$$$

for the diffusive part and

$$$\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$$$

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 \eqref{eq:timedependentweak1} 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.