Galerkin methods
| Slimane Adjerid and Mahboub Baccouch (2010), Scholarpedia, 5(10):10056. | doi:10.4249/scholarpedia.10056 | revision #136217 [link to/cite this article] |
The Galerkin finite element method of lines is one of the most popular and powerful numerical techniques for solving transient partial differential equations of parabolic type. The Galerkin finite element method of lines can be viewed as a separation-of-variables technique combined with a weak finite element formulation to discretize the problem in space. This leads to a stiff system of ordinary differential equations that can be integrated by available off-the-shelf implicit ordinary differential equations (ode) solvers based on implicit time-stepping schemes such as backward differentiation formulas (BDF) or implicit Runge-Kutta (IRK) methods.
Contents |
A one-dimensional parabolic problem
The basic idea of the Galerkin finite element method of lines will be demonstrated on the following one-dimensional linear parabolic partial differential problem
\[\tag{1} \begin{cases} u_t(x,t) - (a(x) u_x(x,t) )_x + b(x) u(x,t) = f(x,t), \ x\in (0,1), \ 0 \le t \le T, \\ u(x,0)=g(x), \ x\in (0,1), \\ u(0,t)=u(1,t)=0, \ 0 \le t \le T, \end{cases} \]
where \(a(x) >0\) and \(b(x) \ge 0 \) on \([0,1]\ .\)
Weak Galerkin formulation
A spatial weak formulation of the model problem is obtained by multiplying the parabolic equation (1) by a test function \(v\) in the Sobolev space \(H_0^1\ ,\) integrating over (0,1) and integrating by parts to show that \(u \in H_0^1 \) satisfies
\[\tag{2} \begin{cases} (v,u_t(.,t)) + A(v,u(.,t)) = (v,f(.,t)), \ \forall \ v \in H_0^1, \ 0 \le t \le T,\\ u(x,0)=g(x), \ x \in (0,1), \end{cases} \]
where
\[ (v,u(.,t)) = \int_0^1 v(x) u(x,t) dx, \ A(v,u(.,t)) = \int_0^1 [a(x) \frac{ {\rm d }v(x)}{{\rm d } x} u_x(x,t) + b(x) v(x)u(x,t)] dx . \]
Semi-discretization
The semi-discrete Galerkin finite formulation consists of subdividing the domain (0,1) into \(n\) elements \(\mathbf{I}_i = (x_i,x_{i+1}), \ i=0,1,\cdots,n \) and constructing a finite-dimensional subspace \(S_p^N\subset H_0^1\ .\) Here \(S_p^N \) consists of piecewise p-degree polynomial functions spanned by a set of basis functions \(\phi_i, \ i=1,2,\cdots,N \ .\) For instance, \( N=n-1, \) for piecewise linear finite elements.
The semi-discrete Galerkin method consists of finding \[ u_N(x,t) = \sum\limits_{i=1}^N c_i(t) \phi_i(x), \] with time-dependent coefficients \(c_i(t) \in \mathbb{R}\) such that
\[ \begin{cases} (v_N,(u_{N})_t(.,t)) + A(v_N,u_N(.,t)) = (v_N,f(.,t)), \ \forall \ v_N \in S_p^N, \ 0 \le t \le T,\\ (v_N, u_N(.,0)=(v_N,g), \ \forall v_N \in S_p^N. \end{cases} \]
Testing against \(v_N = \phi_j, \ j=1,2,\cdots,N\ ,\) leads to the following system of linear ordinary differential equations for \(c_i, \ i=1, \cdots,N\ ,\)
\[\tag{3} \begin{cases} \sum\limits_{i=1}^N \dot{c}_i(\phi_j,\phi_i) + \sum\limits_{i=1}^N c_i A(\phi_j ,\phi_i) = (\phi_j,f(.,t)), \ j=1,2,\cdots,N, \ t \ge 0,\\ \sum\limits_{i=1}^N c_i(0) (\phi_j, \phi_i)=(\phi_j,g), \ j=1,2,\cdots,N. \end{cases} \]
Other optimal approximations of the initial condition include the elliptic projection of
\(g(x)\) onto \(S_p^N\) defined by the equation
\[ A(u_N(.,0),v_N) = A(g,v_N), \ \forall v_N \in S_p^N, \]
and the piecewise polynomial Lagrange interpolant of \(g\) defined by
\[ u_N(x,0) = \sum\limits_{i=1}^N g(x_i) \phi_i(x), \]
where \(\phi_i\) are finite element basis functions of Lagrange type, i.e., \(\phi_i(x_j)=\delta_{ij}\ .\) Here \(\delta_{ij}=1\) if \(i=j\) and 0 otherwise.
Let \( \mathbf{c}(t) = [ c_1(t), c_2(t),\cdots,c_N(t)]^t\ ;\) thus, the linear ordinary differential system (3) may be written in matrix form as
\[\tag{4} \begin{cases} \mathbf{M} \dot{\mathbf{c}} + \mathbf{K} \mathbf{c}(t) = \mathbf{f}(t), \ 0 \le t \le T,\\ \mathbf{M} \mathbf{c}(0) = \mathbf{g} , \end{cases} \]
where the \(N\times N\) mass matrix \(\mathbf M\) and stiffness matrix \(\mathbf K\) defined by \(M_{ij}=M_{ji} =
(\phi_j, \phi_i)\) and \(K_{ij}=K_{ji} = A(\phi_j,\phi_i)\) are symmetric positive definite. The
vectors \( {\mathbf f},\ {\mathbf g} \in \mathbb{R}^N \) are defined by \(f_j(t)= (\phi_j, f(.,t)) \) and
\(g_j=(\phi_j,g)\ .\)
The initial-value problem (4) has a unique solution.
The mass matrix \(\mathbf{M} \) may be replaced by a diagonal matrix using a low-order quadrature to approximate the \(L^2\) inner product, which is easier to invert. This new method is called lumped masses finite element method.
Time stepping
The ordinary differential systems obtained by discretizing parabolic partial differential equations are stiff with a wide disparity in time scales ranging between \(\mathcal{O}(h^{-2})\) and \(\mathcal{O}(1)\ .\) Such systems require the use of implicit and stable time-stepping methods such as backward differentiation formula, implicit Runge Kutta methods (Hairer et al. 1987, Hairer and Wanner 2010) or the popular second-order Crank-Nicolson method to avoid extremely small time steps.
Temporal discretization is obtained by subdividing \([0,T]\) into \(m\) steps of size \(\Delta t =T/m \) and setting \(t_k = k \Delta t, \ k=0,1,\cdots,m \ .\) Then, let \(\mathbf{c}^k\) denote an approximation of \(\mathbf{c}(t_k)\) and write the fully discrete solution at time level \(t_k\) as
\[ u_N^k (x) = \sum\limits_{i=1}^N c_i^k \phi_i(x). \]
The simplest backward differentiation formulas is the backward Euler method which yields
\[\tag{5} \begin{cases} \mathbf{c}^0 = \mathbf{M}^{-1} \mathbf{g}, \\ \mathbf{M} \frac{\mathbf{c}^{k+1} - \mathbf{c}^{k}} {\Delta t } + {\mathbf K} {\mathbf c}^{k+1} = \mathbf{f}(t_{k+1}), \ \ k \ = \ 0,1,\cdots,m-1, \\ \end{cases} \]
while the Crank-Nicolson method leads to
\[\tag{6} \begin{cases} \mathbf{c}^0 = \mathbf{M}^{-1} \mathbf{g}, \\ \left ( \frac{\mathbf{M}}{\Delta t} + \frac{\mathbf{K}}{2} \right )\mathbf{c}^{k+1} = \left ( \frac{\mathbf{M}}{\Delta t} - \frac{\mathbf{K}}{2} \right ) \mathbf{c}^k + \mathbf{f}(\frac{t_{k+1}+t_{k}}{2}), \ \ k \ = \ 0,1,\cdots,m-1. \end{cases} \]
Both methods require the solution of a linear algebraic system at each step to compute \(\mathbf{c}^{k+1}\ .\)
The discontinuous Galerkin method in time is stable and equivalent to implicit Radau Runge-Kutta methods (Karakashian 1998). Thus, it may be applied as a time-stepping method to solve the ode system (4).
An example
Here we use a uniform mesh with \(h=1/n\ ,\) \(x_i =i h, \ i=0,1,\cdots,n\) and construct a piecewise linear basis \(\phi_i, \ i=1,\cdots,n-1,\) for \(S_1^{n-1}\ .\) A basis function \(\phi_i\) is defined below with its graph shown in Figure 1.
\[ \phi_i(x) = \begin{cases} \frac{ x-x_{i-1}}{h}, & x_{i-1} \le x \le x_i, \\ \frac{ x_{i+1}-x}{h}, & x_{i} \le x \le x_{i+1}, \\ 0, & elsewhere. \end{cases} \]
If \( a=1, b=0\) and \(f =0\ ,\)
\[ \mathbf{M} = \frac{h}{6}\begin{bmatrix} 4 & 1 & 0 & \cdots & 0 \\ 1 & 4 & 1 & \ddots & 0 \\ 0 & 1 & 4 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots& 0 & 1 & 4 \end{bmatrix}, \quad \mathbf{K} = \frac{1}{h} \begin{bmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \ddots & 0 \\ 0 & -1 & 2 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -1 \\ 0 & \cdots & 0 & -1& 2 \end{bmatrix}, \quad \mathbf{F} = \mathbf{0}, \ \mathbf{c}^0 = \begin{bmatrix} g(x_1) \\ g(x_2) \\ g(x_3) \\ \vdots \\ g(x_{n-1}) \end{bmatrix}. \]
Instead of \(\mathbf M\ ,\) one may use the diagonal lumped mass matrix, \(\mathbf M=h\mathbf I\ ,\) obtained using the 2-point Gauss-Lobatto quadrature with \(\xi=[-1,1]\) and weights \(w=[1,1]\ .\)
Multidimensional problems
A Galerkin finite element method for multi-dimensional parabolic problems is derived in a similar manner, as illustrated on the following model problem:
\[\tag{7} \begin{cases} u_t(\mathbf{x},t)=\Delta u (\mathbf{x},t)+f({\mathbf{x}},t), \quad ({\mathbf x},t) \in \Omega\times[0,T],\\ u(\mathbf{x},t) = 0, \quad (\mathbf{x},t) \in \partial \Omega \times [0,T],\\ u(\mathbf{x},0) = g(\mathbf{x}) , \qquad \mathbf{x} \in \Omega, \end{cases} \]
where \(\Omega \subset \mathbb{R}^d \ .\)
The domain \(\Omega \) is partitioned into a regular mesh consisting of triangular or rectangular elements (in two dimensions) and tetrahedral or hexahedral elements in three dimensions. We apply the standard procedure to construct a finite-dimensional subspace \(S_p^N\) of \(H_0^1\) spanned by a piecewise p-degree polynomial finite element basis \(\phi_i(\mathbf{x}), \ i=1,2,\cdots,N,\) developed for elliptic problems; see (Axelson and Barker 2006, Baker et al. 1981).
The semi-discrete Galerkin finite element problem consists of finding
\[ u_N(\mathbf{x},t) = \sum\limits_{i=1}^N c_i(t) \phi_i(\mathbf{x}), \]
such that
\[ \begin{cases} (v_N,(u_{N})_t(.,t)) + A(v_N,u_N(.,t)) = (v_N,f(.,t)), \ \forall \ v_N\in S_p^N, \ t \ge 0,\\ (v_N, u_N(.,0)=(v_N,g), \ \forall \ v_N \in S_p^N, \end{cases} \]
where
\[ (v,u(.,t))= \int_{\Omega} v(\mathbf{x}) u(\mathbf{x},t) d \mathbf{x}, \ A(v,u(.,t)) = \int_{\Omega} [ \nabla v(\mathbf{x})]^t \nabla u (\mathbf{x},t) d\mathbf{x}. \]
Testing against \(v_N = \phi_j, \ j=1,2,\cdots,N\) yields a system of differential equations of the form (4).
Again, an implicit time-stepping method should be used to integrate the ordinary differential system in time.
A priori error estimates
If the exact solution has enough smoothness so that the terms appearing in the right-hand side of the error bounds given below are finite and if \(h_k\) denotes the diameter of element \(\Delta_k\) and \(h= \max\limits_k h_k \ ,\) then the following a priori error estimates on a family of regular meshes hold.
Semi-discrete errors
There are well known optimal a priori error estimates for the semi-discrete solution in the \(L^2\) norm as in (Thomee 2006).
The error between the true and the semi-discrete solutions in the \(L^2\) norm is bounded as:
\[ || u(.,t)-u_N(.,t)|| \leq || g -u_N(.,0)||+ C h^{p+1} ( || g ||_{p+1} + \int_0^t ||u_t(.,\tau)||_{p+1} d\tau ), \ 0 \le t \le T, \]
where \( || \cdot || \) and \(|| \cdot ||_{s}\ ,\) respectively, denote \(L^2\) and \(H^s\) Sobolev norms. Here and in all the estimates below, \(C \) is a generic positive constant independent of \(h\ .\)
The error between the gradients of the true and the semi-discrete solutions in the \(L^2\) norm is bounded as:
\[ ||\nabla u(.,t)-\nabla u_N(.,t)|| \leq C || \nabla u_N(.,0)-\nabla g || \] \[ + C h^{p} \left ( || g ||_{p+1} + ||u(.,t)||_{p+1} + [\int_0^t ||u_t(.,\tau)||_p^2 d \tau ]^{1/2} \right ), \ 0 \leq t \leq T. \]
For smooth solutions and optimal approximations of initial conditions the semi-discrete error is bounded by \[ || (u -u_h) (.,t)|| + h || (\nabla u - \nabla u_h)(.,t)|| \le C(u) h^{p+1}. \]
If in addition the mesh is quasi-uniform, the error between the true and the semi-discrete solutions in the \(L^\infty\) norm for linear elements, \(p=1\ ,\) is almost optimal and is given by
\[ \max\limits_{t \in [0, T]} || u(.,t - u_N(.,t) ||_{L_\infty} \leq C (1 + log(\frac{1}{h})) h^2 \max\limits_{t \in [0,T]} || u(.,t)||_{W^2_\infty}. \]
Fully discrete errors
Using the backward Euler time-stepping scheme, the error between the true solution and the fully discrete solution in the \(L^2\) norm can be bounded as
\[ || u(.,t_k) - u_N^k || \le C(u) (h^{p+1} + \Delta t), \ C(u) >0. \]
The \(L^2\) error between the true solution and the fully discrete solution with linear finite elements, \(p=1\ ,\) and the Crank-Nicolson method is given as
\[ || u(.,t_k) - u_N^k || \le C(u) (h^{p+1} + \Delta t ^2), \ C(u) >0. \]
References
- O. Axelsson and V.B. Barker., Finite Element Solution of Boundary Value Problems, SIAM, Philadelphia, (2001).
- E.B. Baker, G.F. Carey and J.T. Oden, Finite Elements, An Introduction, Volume I, Prentice Hall, Engelwood Cliffs, New Jersey, (1981).
- E. Hairer, S. P. Norsett and G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, Springer Verlag, Berlin, (1987).
- E. Hairer and G. Wanner, Solving Ordinary Differential Equations I, Stiff and Differential-Algebraic Problems, Springer Verlag, Berlin, (2010).
- O. Karashian and C. Makridakis, A space-time finite element method for the nonlinear Schrodinger equation: The discontinuous Galerkin method, Mathematics of Computation, 67, pp. 479-499. (1998).
- V. Thomee, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, Berlin, (2006).
Internal references
- Bill Gear (2007) Backward differentiation formulas. Scholarpedia, 2(8):3162.
- J. Tinsley Oden (2010) Finite element method. Scholarpedia, 5(5):9836.
- Andrei D. Polyanin, William E. Schiesser, Alexei I. Zhurov (2008) Partial differential equation. Scholarpedia, 3(10):4605.
- John Butcher (2007) Runge-Kutta methods. Scholarpedia, 2(9):3147.
- Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.
- Lawrence F. Shampine and Skip Thompson (2007) Stiff systems. Scholarpedia, 2(3):2855.
See also
(pending)
Prof. Slimane Adjerid accepted the invitation on 21 September 2009 (self-imposed deadline: 21 March 2010).
