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∈Vu \in V such that

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

for all v∈Vv \in V and t∈[0,∞)t \in [0,\infty).

We can convert (1) into a system of equations

M(t)u˙+A(t)u=b(t)\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

MI ⁣J(t)=m(NJ,NI;t)AI ⁣J(t)=a(NJ,NI;t)bI(t)=b(NI;t).\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 u˙≈[un+1−un]/Δt\dot{u} \approx [u_{n+1}-u_n]/{\Delta t}, linearity allows us to write

m(u˙,v;t)≈1Δt[m(un+1,v;tn+1)−m(un,v;tn)]\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

M(t)u˙≈1Δt[Mn+1un+1−Mnun]\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, un+1:=u(x,tn+1)u_{n+1}:= u(x, t_{n+1}), Mn+1=M(tn+1)\BM_{n+1} = \BM(t_{n+1}) and vice versa for unu_n and Mn\BM_n.

Explicit Euler Scheme

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

1Δt[m(un+1,v;tn+1)−m(un,v;tn)]+a(un,v;tn)=b(v;tn)∀v∈V .\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

1Δt[Mn+1un+1−Mnun]+Anun=bn\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

un+1=Mn+1−1[Mnun+Δt(bn−Anun)]\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 mm is time-independent, that is m(u˙,v;t)=m(u˙,v)m(\dot{u}, v;t) = m(\dot{u}, v), we have

un+1=un+Δt M−1(bn−Anun)\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 tn+1t_{n+1}

1Δt[m(un+1,v;tn+1)−m(un,v;tn)]+a(un,v;tn+1)=b(v;tn+1)∀v∈V .\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

1Δt[Mn+1un+1−Mnun]+An+1un+1=bn+1\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

un+1=[Mn+1+ΔtAn+1]−1[Mnun+Δt bn+1]\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 mm is time-independent, one can just substitute M=Mn+1=Mn\BM=\BM_{n+1}=\BM_n.

Example: Reaction-Advection-Diffusion Equation

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

∂u∂t=∇⋅(D∇u)−∇⋅(cu)+ru+finΩ×Iu=uˉon∂Ω×Iu=u0inΩ,t=0\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∈I=[0,∞)t\in I = [0,\infty),

  • D\BD is a second-order tensor describing the diffusivity of uu,
  • c\Bc is a vector describing the velocity of advection,
  • rr is a scalar describing the rate of reaction,
  • and ff is a source term for uu.

The weak formulation is then

Find u∈Vu \in V such that

∫Ωu˙v dv=∫Ω[∇⋅(D∇u)−∇⋅(cu)+ru+f]v dv\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∈Vv \in V and t∈It \in I.

We have the following integration by parts relationships:

∫Ω∇⋅(D∇u)v dv=∫Ω∇⋅(vD∇u) dv−∫Ω(D∇u)⋅∇v dv\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

∫Ω∇⋅(cu)v dv=∫Ω∇⋅(cuv) dv−∫Ωuc⋅∇v dv\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=0v=0 on the boundary. Then our variational formulation is of the form (1) where

m(u˙,v)=∫Ωu˙v dva(u,v)=∫Ω(D∇u)⋅∇v dv−∫Ωuc⋅∇v dv−∫Ωruv dvb(v)=∫Ωfv dv\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

MI ⁣J=∫ΩNJNI dvAI ⁣J=∫Ω(DBJ)⋅BI dv−∫ΩNJc⋅BI dv−∫ΩrNJNI dvbI=∫ΩfNI dv\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 M\BM is constant through time.