Ngô Quốc Anh

February 24, 2010

FEM for heat equations in 2D

Filed under: Giải tích 9 (MA5265) — Ngô Quốc Anh @ 19:01

In this topic, we have introduced the finite element method via the Poisson equations in 1D. Note that the FEM has some advantages compared to the FDM, for example, the mesh is not necessarily uniform.  Today, we discuss the question how to apply this FE method to some evolution equations. What we are going to do is to study the heat equations in 2D.

In the literature, the simplest model of the heat equations in 2D over a domain \Omega is given as follows

\displaystyle \frac{\partial u}{\partial t} - \frac{\partial^2 u}{\partial x_1^2} - \frac{\partial^2 u}{\partial x_2^2} = f(t,x_1,x_2)

where f is a given function. For simplicity, we write u_t - \Delta u = f(t,x). Unlike the heat equations in 1D, we cannot illustrate the solution u in just one picture some several t. The point is u has 3 variables, and thus to achieve a full picture of u, we need to draw u in 4D, that’s impossible.

Due to the presence of t, normally we need to impose an initial condition, i.e., the solution u at t=0. The way to solve this equation numerically is to use, for example, the Backward Euler Method. We denote by dt the difference in time, if we have already known the solution u at t=t_0, then we can solve u at t=t_0+dt. The Backward Euler Method applied to this model equation gives

\displaystyle\frac{{u({t_0} + dt,x) - u({t_0},x)}}{{dt}} - \Delta u({t_0} + dt,x) = f({t_0} + dt,x).

Usually, for the sake of clarity by n we mean we are working on the n-th step, therefore, at n-step, the time t is given by t_n=ndt, the solution u is then denoted by u^n(x) which is actually u(t_n, x). The time-discretization equation is then given by

\displaystyle\frac{{{u^n}(x) - {u^{n - 1}}(x)}}{{dt}} - \Delta {u^n}(x) = f({t_n},x).

We assume T is a triangulation of \Omega, that means, we devide the domain \Omega into M small special triangles. Vertices of triangles, usually called nodes, are those points we want to approximate u. If we denote by V_h the finite-dimensional subspace of H^1(\Omega) with a basis \{\varphi_1,...,\varphi_M\} such that for each u_h \in V_h, the restriction of u_h into T is piecewise linear in the sense that u_h is linear on each triangle of T.

The above picture shows us a case when \Omega is a ball. The idea of FEM is to approximate solution u by some u_h \in V_h such that u(x)=u_h(x) for every vertices x. Let talk about the basis \{\varphi_1,...,\varphi_M\}. Functions \varphi_i at the i-node are chosen so that \varphi({\rm node}_i)=1 and \varphi({\rm node}_j)=0 for every j \ne i.

We are now interested in finding weak solutions to the heat equations. In this situation, one tries to find a function u_h \in V_h which satisfies the following equation

\displaystyle\frac{1}{{dt}}\int_\Omega {u_h^n(x)v(x)dx} - \frac{1}{{dt}}\int_\Omega {u_h^{n - 1}(x)v(x)dx} + \int_\Omega {\nabla u_h^n(x)\nabla v(x)dx} = \int_\Omega {f({t_n},x)v(x)dx}

for every v\in V_h.


\displaystyle u_h^i(x) = \sum\limits_{k = 1}^M {{\xi ^i}{\varphi _k}(x)}

we then have

\displaystyle\begin{gathered} \sum\limits_{k = 1}^M {{\xi ^n}\left( {\int_\Omega {{\varphi _k}(x)v(x)dx} + dt\int_\Omega {\nabla {\varphi _k}(x)\nabla v(x)dx} } \right)} \hfill \\ \qquad\qquad= \sum\limits_{k = 1}^M {{\xi ^{n - 1}}\int_\Omega {{\varphi _k}(x)v(x)dx} } + dt\int_\Omega {f({t_n},x)v(x)dx} \hfill \\ \end{gathered}.

If we apply the above equation for v(x) being \varphi_j(x), we then have

\displaystyle\begin{gathered} \sum\limits_{k = 1}^M {{\xi ^n}\left( {\int_\Omega {{\varphi _k}(x){\varphi _j}(x)dx} + dt\int_\Omega {\nabla {\varphi _k}(x)\nabla {\varphi _j}(x)dx} } \right)} \hfill \\ \qquad\qquad= \sum\limits_{k = 1}^M {{\xi ^{n - 1}}\int_\Omega {{\varphi _k}(x){\varphi _j}(x)dx} } + dt\int_\Omega {f({t_n},x){\varphi _j}(x)dx} , \qquad j = \overline {1,M} . \hfill \\ \end{gathered}.

In term of linear algebra, the above equations can be transfered to a linear system (A+dtB) \xi^n = B\xi^{n-1}+b where

\displaystyle B = \left( {\begin{array}{*{20}{c}} {\int_\Omega {{\varphi _1}(x){\varphi _1}(x)dx} } & \cdots & \cdots & {\int_\Omega {{\varphi _1}(x){\varphi _M}(x)dx} } \\ \vdots & \ddots & {} & \vdots \\ \vdots & {} & \ddots & \vdots \\ {\int_\Omega {{\varphi _M}(x){\varphi _1}(x)dx} } & \cdots & \cdots & {\int_\Omega {{\varphi _M}(x){\varphi _M}(x)dx} } \\ \end{array} } \right)


\displaystyle A = \left( {\begin{array}{*{20}{c}} {\int_\Omega  {{\nabla\varphi _1}(x){\nabla\varphi _1}(x)dx} } & \cdots & \cdots &  {\int_\Omega {{\nabla\varphi _1}(x){\nabla\varphi _M}(x)dx} } \\ \vdots & \ddots  & {} & \vdots \\ \vdots & {} & \ddots & \vdots \\  {\int_\Omega {{\nabla\varphi _M}(x){\nabla\varphi _1}(x)dx} } & \cdots &  \cdots & {\int_\Omega {{\nabla\varphi _M}(x){\nabla\varphi _M}(x)dx} } \\  \end{array} } \right)


\displaystyle b = dt\left( {\begin{array}{*{20}{c}} {\int_\Omega {f({t_n},x){\varphi _1}(x)dx} } \\ \vdots \\ \vdots \\ {\int_\Omega {f({t_n},x){\varphi _M}(x)dx} } \\ \end{array} } \right).

The rest of our task is to solve this system of linear equations.

Source: and a book due to Claes Johnson.


Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

Create a free website or blog at

%d bloggers like this: