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

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

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

and $\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: http://en.wikipedia.org/wiki/Finite_element_method and a book due to Claes Johnson.