Solving Parabolic Problem with FEM

The method of finite elements can be applied also to initial value problems, such as the parabolic problem. For an alternative numerical method, the finite difference method, the main theory is well known: A formally consistent method is convergent if and only if it is stable. This statement holds also for FEM, and accoring to Strang, the only difference is that stability and consistency are often easier to verify for finite elements than for finite differences.

In this presentation some solution methods are applied to a simple heat equation. First the characteristics of the problem are considered, then the spatial discretization using FEM is conducted. Next the classical methods and the discontinuous Galerkin approach for the time discretization are introduced.

This text is a brief presentation of the ideas brought up by Gilbert Strang and George Fix with their book "An analysis of the finite element method" (1973) and Claes Johnson with his book "Numerical solution of partial differential equations by the finite element methods" (1987).

Background

Let us first consider a one-dimensional problem of modelling heat conduction in a bar:

\frac{∂u}{∂t} − \frac{∂^{2}u}{∂x^{2}} = f(x, t),  0 < x < π, t > 0  (1)

with the following boundary and initial condition:

u(0, t) = u(π, t) = 0, t > 0;  u(x, 0) = u^{0}(x), 0 < x < π.  (2)

If the source-term f is taken to be zero, the solution of the problem is given by

u(x, t) = ∑_{j=1}^{∞} u_{j}^{0} e^{−j^{2}t} sin(jx)  (3)

where

u_{j}^{0} = \sqrt{2/π} ∫_{0}^{π} u^{0}(x) sin(jx) dx,  j = 1, 2, ….  (4)

The terms u_{j}^{0} are actually the Fourier coefficients of the initial data u^{0}. Thus u(x, t) is a linear combination of the sine waves sin(jx) with frequencies j and amplitudes u_{j}^{0} exp(−j^{2}t). For a reasonably large j^{2}t the term exp(−j^{2}t) is very small, and therefore high frequency components will be dampened quickly. Hence it can be said that each component sin(jx) lives on a time scale of order O(j^{−2}).

Due to the dampening, the solution u(x, t) will become smoother as t increases. In general, for small t the solution will not be smooth, because the sizes of the derivatives of u with respect to x or t depend on how quickly the Fourier coefficients decay. This initial transient requires the use of smaller mesh size both in time and space. It should be noted that transients may occur also if for example the right hand side f or the boundary conditions vary sharply in time.

Spatial Discretization

Discretizing the problem in space only yields an initial value problem for a system of ordinary differential equations. To begin with, let us write the weak formulation of the problem:

∫_{0}^{π}(u_{t} − u_{xx} − f)v dx = 0.  (5)

For greater symmetry between the trial function and test function, the term −u_{xx}v is integrated by parts, yielding

∫_{0}^{π} u_{t}v + u_{x}v_{x} − fv dx = 0.  (6)

Now given an N-dimensional subspace S^{h} of H_{E}^{1}, the Galerkin principle is to find a function u^{h}(x, t) with the following property: At each t > 0, u^{h} lies in S^{h} and satisfies

∫_{0}^{π} u_{t}^{h}v^{h} + u_{x}^{h}v_{x}^{h} − fv^{h} dx = 0 ∀ v^{h} ∈ S^{h} .  (7)

The Galerkin formulation is discrete in space variables but continuous in time variable. Choosing a basis φ_{1}, …, φ_{N} for the trial space S^{h} the unknown solution can be expanded as

u^{h}(x, t) = ∑_{1}^{N} Q_{j}(t) φ_{j}(x),  (8)

where the optimal weights Q_{j} are determined by the Galerkin principle:

∑_{j} ∫_{0}^{π}(\frac{∂Q_{j}}{∂t} φ_{j}φ_{k} + Q_{j} \frac{∂φ_{j}}{∂x} \frac{∂φ_{k}}{∂x} − fφ_{k}) dx = 0,  k = 1, …, N.  (9)

The result is a system of N ordinary differential equations for the N unknowns Q_{1}(t), …, Q_{N}(t); the boundary conditions are already incorporated. An efficient way to account for the initial condition is to use the interpolate of u^{0}: u_{0}^{h} = (u^{0})_{I} .

Writing the Galerkin equations in matrix form, we get:

MQ′ + KQ = F (t) ,  (10)

where M is the mass matrix, K is the stiffness matrix — the same ones that appear in the steady-state problem. The right hand side is F_{k} = ∫ f(x, t) φ_{k}(x) dx .

Classical methods

The obtained system of ordinary differential equations is usually stiff: the solution (Eq. 3) contains components with vastly different time scales. Formally this can be seen also from the considerable variation in the size of the eigenvalues of the prolem.

For the time-discretization of the system there are three basic methods: explicit Euler’s method, implicit Euler’s method (backward Euler) and trapezoidal rule (Crank-Nicolson). In order to an explicit method to track the characteristics of a stiff system, especially the initial transient, the time step k_{n} has to be very small, or more precisely

k_{n} ≤ Ch^{2}  (11)

where h is the mesh size. Another way to put this is to say that k_{n} has to be chosen so small that the fastest time scale is resolved. The implicit methods are more efficient, because after the initial transient phase the time step can be increased. They are also unconditionally stable, unlike the explicit methods.

Replacing the time derivatives with approximations in accordance with the Crank-Nicolson scheme yields

M \frac{Q^{n+1} − Q^{n}}{Δt} + K \frac{Q^{n+1} + Q^{n}}{2} = \frac{f^{n+1} + f^{n}}{2} ,  (12)

which can be rewritten as

Q^{n+1} = (\frac{M + KΔt}{2})^{−1} (\frac{M − KΔt}{2} Q^{n} + \frac{Δt}{2}(f^{n+1} + f^{n})) .  (13)

It automatically filters out the useless high modes if the time step is chosen properly.

The implicit Euler’s method leads to quite similar result, but it is even better than Crank-Nicolson: The latter is A-stable when the former is L-stable. One often used trick is to calculate a couple of steps with the non-oscillating Implicit Euler and then continue with the more accurate Crank-Nicolson scheme.

Discontinuous Galerkin method

A question may arise that way FEM is not used for the temporal discretization. It can be used, but for example Strang argues, that in the time variable there are no geometrical difficulties to overcome by finite elements. He continues that the straight forward application of the Galerkin principle may couple all the time levels and distroy the crucial property of propagation forward in time. Therefore he sees no reason the give up the flexibility of finite difference methods and concludes, that the attempts done have not been successfull.

However, Claes Johnson considers a more recent method that is based on finite element formulation in time with piecewise polynomials of degree q. He argues that the advantage of this method is that the variable coefficients and non-zero right hand sides present no complications in principle. Further the variational formulation is very useful in the analysis of the time discretization error. He admits that the possibilities offered by the discontinuous Galerkin method have been dicovered only recently and they are still under exploration.

Johnson presents the discontinuous Galerkin method as follows. After the spatial discretization we have:

(u̇_{h}(t), v) + a(u_{h}(t), v) = (f(t), v) ∀v ∈ V, t ∈ I; u(0) = u^{0}  (14)

where for our simple case a = 1. Now in order to formulate this further, a new space for a given non-negative integer q is introduced:

W_{hk} = {v :→ V_{h} : v|_{I_{n}} ∈ P_{q}(I_{n}), n = 1, …, N},  (15)

where

P_{q}(I_{n}) = {v :→ V_{h} : v(t) = ∑_{i=0}^{q} v_{i} t^{i} with w_{i} ∈ V_{h}}.  (16)

In other words, W_{hk} is the space of functions on I with values in V_{h} that on each time interval I_{n} vary as polynomials of degree at most q. It should be noted that the functions v in W_{hk} may be discontinuous in time at the discrete time levels t_{n}. For this reason, let us use the notation

v_{+}^{n} = lim_{s→0^{+}} v(t_{n} + s), v_{−}^{n} = lim_{s→0^{−}} v(t_{n} + s), [v^{n}] = v_{+}^{n} − v_{−}^{n} ,  (17)

where [v^{n}] is the jump of v at t_{n}.

Now the discontinuous Galerkin method can be expressed as follows: Find U ∈ W_{hk} such that

A(U, v) = L(v),  ∀v ∈ W_{hk}  (18)

where

A(w, v) = ∑_{n=1}^{N} ∫_{I_{n}}(ẇ, v) + a(w, v) dt + ∑_{n=2}^{N}([w^{n−1}], v_{+}^{n−1}) + (w_{+}^{0}, v_{+}^{0})  (19)
L(v) = ∫_{I}(f, v)dt + (u^{0}, v_{+}^{0}) .  (20)

Since v ∈ W_{hk} varies independently on each subinterval I_{n}, there is an alternative formulation: For n = 1, …, N, given U_{−}^{n−1}, find U ≡ U|_{I_{n}} ∈ P_{q}(I_{n}) such that

∫_{I_{n}} ((U̇, v) + a(U, v)) dt + (U_{+}^{n−1}, v_{+}^{n−1}) = ∫_{I_{n}}(f, v) dt + (U_{−}^{n−1}, v_{+}^{n−1})  (21)

∀v ∈ P_{q}(I_{n}) and where U_{−}^{0} = u^{0}.

In simple cases the discontinuous Galerkin method gives schemes which coincide with classical method based on so called sub-diagonal Padé approximations. For example, choosing q = 0, and using the notation U^{n} ≡ U_{−}^{n} ≡ U_{+}^{n−1}, the formulation above reduces to finding U^{n} ∈ V_{h} such that

(U^{n} − U^{n−1}, v) + k_{n}a(U^{n}, v) = ∫_{I_{n}}(f, v) dt ∀v ∈ V_{h}, n = 1, …, N ,  (22)

where U_{−}^{n} = u^{0}. This is a modification of Implicit Euler’s method, where the right hand side involves an average of f over I_{n} rather than the value of f at t_{n}.