FEM4EllipticPDE/dox/chapters/derivation.tex

\chapter{Derivation}

\section{Problem Statement}
FEM4EllipticPDE is a Lemma module implementing a Galerkin finite element solution to a family of commonly encountered 
elliptic  problems taking the form
%Considering the uniformly elliptic Dirichlet BVP: 
\begin{eqnarray}
	\label{eq:scheme}
	-\nabla \cdot \left( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) \right) = g(\mathbf{r}) & \mathop{\forall}_{u \in \Omega} \\
	u = 0                                            & \mathop{\forall}_{u \in \partial \Omega} 
\end{eqnarray}

Where $0 < \sigma_{min} \leq \sigma(\mathbf{r}) \leq \sigma_{max} \ll \infty $  and $\Omega \in \mathbb{R}^3$. When $\sigma \equiv 1$, this 
reduces to Poisson's equation $\left( \nabla^2 u  = g \right)$. These types of problem arises in several areas of geophysics: gravitational or magnetic potentials 
obey Poisson's equation, while the more general form can be used to solve electrostatic (DC and SP) problems. 

\begin{enumerate}
	\item The Galerkin finite element method is defined by reposing the global problem in terms of numerous local ones through the use of 
	appropriate test functions. 
	\begin{itemize}

		\item Test functions $v$ must be from a suitable subspace. 
		In this case, due to the Dirichlet boundary conditions,
		\[
			v  \in H^a_b(\Omega)
		\]
		meet our requirements. $H^a_b$ represents the subspace of weakly differentiable functions that are zero valued at $a$, and $b$
		and $\textgoth{V}$ represents the actual functions. % Hilbert space.

		The 3D variational problem may be written

%		\[
%				-\nabla \cdot \left( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) \right) = g(\mathbf{r}) & \mathop{\forall}_{u \in \Omega} \\
%	u = 0                                            & \mathop{\forall}_{u \in \partial \Omega} 
%
		%\] 
		\[
		  \int	-v \nabla \cdot \sigma  \nabla u = v g.
		\]
		This can be simplified 
		\[
			\int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx = v g
		\]
		The proof of which follows below. 
		%\item The variational form may be found by multiplying the elliptic PDE by any arbitrary
		%\[ 
		%	v \in H_a^b(\Omega) 
		%\]
		
		%In 1D ($\Omega \in \mathbb{R}^1$) \autoref{eq:scheme} reduces to:
		%\begin{equation} \label{eq:oned}
		%	- \DD{x} \sigma \DD{x} u = g
		%\end{equation}
		%Using the variational (weak) formulation for any arbitrary $v \in \textgoth{V} 
	
		\item In any arbitrary dimension $x$, if the boundaries in that dimension of $\partial \Omega = [a,b]$ the variational problem using 
		test functions can be constructed 
		\begin{equation} \label{eq:oned}
			-\int_a^b v(x) \DD{x} \sigma(x) \DD{x} u(x)  dx = \int_a^b v(x) g(x) dx ,  
		\end{equation}

				
		The variation problem can be reduced to a simplified version, if appropriate test functions we may instead write the 
		variational problem as 
		\[
			\int_a^b  \sigma \frac{du}{dx} \frac{dv}{dx} dx = \int_a^b v(x) g(x) dx
		\]
		

		In one dimension this resolves to
		\[
			%	-\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx = 
			\int_a^b  \sigma \frac{du}{dx} \frac{dv}{dx} dx = \int_a^b v(x) g(x) dx
		\]

		
				

		\begin{corollary}
		\[
		\frac{d}{dx} \left[ v \sigma \frac{du}{dx}  \right] = v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] +
		 \sigma \frac{du}{dx} \frac{dv}{dx}  
		\]	
		\end{corollary}

		

		\begin{proof}
		and the integration by parts on the left hand side yields
		\begin{eqnarray*}
		\int_a^b v \frac{d}{dx} \left[ v \sigma \frac{du}{dx} \right] dx &=&  \int_a^b v \frac{d}{dx} 
		\left[ \sigma \frac{du}{dx} \right] dx + 
		\int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx  \\
		\Rightarrow & & \\
		\left[ v \sigma \frac{du}{dx} \right]_a^b &=& \int_a^b v \frac{d}{dx} \left[ 
		\sigma \frac{du}{dx} \right]  + \int_a^b \sigma \frac{dv}{dx}\frac{du}{dx} dx \\
		\Rightarrow & & \\\
			-\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx &=& \left[ v \sigma 
			\frac{du}{dx} \right]_a^b + \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx
		\end{eqnarray*}
		\[
				-\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx = \int_a^b 
		\sigma \frac{du}{dx} \frac{dv}{dx} dx
		\]
		\end{proof}

		Because $v\in H^1_0$ the solution vanishes at the boundaries $\left[ v \sigma \frac{du}{dx} \right]_a^b =0 $.
		Therefore the scheme reduces to:

% 		Which can be rewritten as
% 		\[
% 			a(u,v) = <g,v>
% 		\] 
		
		The Galerkin FEM method uses triangle (hat) functions for $v$, which are renamed $\phi$ as 
		they are specific. I will occasionally use these interchangeably.  

		Through induction it can be shown then that the 3D problem can similarity be posed 
		So that Equation \ref{eq:scheme} may be rewritten
		\[
			-v \nabla \cdot \sigma  \nabla u = g v.
		\]

		\FloatBarrier
	
		\item \FloatBarrier Define mesh. Since $\Omega = (0,\pi)$ is an uncountable subset of $\mathbb{R}^1$ a 
		countable subset of $\Omega$ must be defined. A mesh is defined over the interval $(0,\pi)$.
		This mesh is defined and shown in Figure \ref{fig:mesh}.	

		\begin{figure}[ht]
		\begin{center}
		\input{chapters/mesh}
		\caption{A 1D FEM mesh is defined. $N$ nodes ($n$) are defined over the interval. Between the 
		nodes an element is defined. The mesh is general and the spacing between nodes is defined by 
		$N-1$ discritisation parameters $h$.}
		\label{fig:mesh}
		\end{center}
		\end{figure}

		\item Define a finite dimensional subspace based on the mesh and a particular test function.  
		A subspace  $\textgoth{V}_h \subset \textgoth{V}$ is sought such that $ \textgoth{V}_h$ is finite 
		dimensional. 
		Particular linearity independent test functions $\phi$ are defined on $\textgoth{V}_h$ such that 
		$\textgoth{V}_h = \mathrm{span}(\phi_i)$. The test functions form a basis of this space. 
		If the weak formulation of the derivative is used these
		test functions may be constructed such as they are only piecewise differentiable. We may define test 
		functions as Kronecker delta $\delta$ functions taking values of unity at a particular node of the mesh. 		
		The linear spline Galerkin method of FEM interpolates these $\delta$ functions. They are zero valued at the 
		neighbouring nodes. In 1D this interpolation yields hat functions.
		Therefore the test functions  $\phi_i$  are defined as:

		\begin{equation}
		\phi_i(x) = \left\{ \begin{array}{ccl} 
		0 & \textrm{if}  & x\in[\alpha, x_{i-1}]  \\
		\frac{x - x_{i-1}}{h_i} & \textrm{if}  & x\in[x_{i-1}, x_i]  \\
		\frac{x_{i+1}-x}{h_{i+1}} & \textrm{if}  & x\in[x_i, x_{i+1}]  \\
		0 & \textrm{if}  & x\in[x_{i+1}, \beta ]  \\
		\end{array} \right.
		\end{equation} 

		Define the variational form of the solution solution 

		\begin{equation} \label{eq:sol1}
		a(u, v) = <g,v> 
		\end{equation}

		A is a symmetric bilinear form such that $a(u,v) = a(v,u)   \forall v \in \textgoth{V}_h$.
		
	\end{itemize}

	In 3D things are a little different. For a given element (tetrahedra) comprised of 4 points we construct the 
	location matrix
	\[ \mathbf{C}  = 
		\left[ \begin{matrix}
			1 & x_1 & y_1 & z_1  \\
			1 & x_2 & y_2 & z_2  \\
			1 & x_3 & y_3 & z_3  \\
			1 & x_4 & y_4 & z_4
		       \end{matrix} \right]
	\]

	We can then fill the stiffness matrix
	\[
		K_{ij} = \sum_{i=0}^{N} \int_{Tet_k} = c \nabla \phi_i \cdot \phi_j
	\]

	\item The solution may be written in matrix(stiffness matrix) vector (load vector) notation. 

	We approximate the solution (Equation \ref{eq:sol1}) using the trial and test functions $\phi$ discussed 
	previously. Let 
	\begin{eqnarray*}
		\sigma(x) &=& x + 1 \\
		u(x) &=& \sin(x) 
	\end{eqnarray*}

	At each node $u_h \in \textgoth{V}_h : a(u_h, v) = <g,v>$ the stiffness matrix $A$
	\begin{eqnarray*}
		A_{ij} &=& a(\phi_i, \phi_j) = \int_0^\pi \sigma(x) \phi'_i \phi'_j \\
		       &=&  \int_0^\pi (x+1) \phi'_i \phi'_j 
	\end{eqnarray*}
	
	%//Since $\phi$ is a triangle, $\phi'$ becomes a step function.

	We have for $i= 1,\cdots , n$
	\[
		\phi'(x) = \left\{ \begin{array}{ccl}
			0 & \textrm{if}  & x\in[\alpha, x_{i-1}]  \\
			\frac{1}{h_i} & \textrm{if}  & x\in[x_{i-1}, x_i]  \\
			-\frac{1}{h_{i+1}} & \textrm{if}  & x\in[x_i, x_{i+1}]  \\
			0 & \textrm{if}  & x\in[x_{i+1}, \beta ]  \\
		\end{array} \right.
	\]

	The product $\phi_i'(x) \phi_j(x)$ vanishes when $|i - j| > 1$. So that in a 1D case, A becomes tri-diagonal. 
	Following the general equation:
	\begin{eqnarray*}
		A_{ij} &=&  \int_{i-1}^{i+1} (x+1) \phi'_i(x) \phi'_j(x) \\ 
	\end{eqnarray*}

	The load vector $g$ is given:
	\[
		[g]_i = <g, \phi_i> = \int_0^\pi g \phi_i
	\]

	Using $ g = [g]$, the linear system becomes:
	\[
		Au = g
	\]

	\item For this case the entries of the stiffness matrix  and load vector are
	
	For the three non-zero cases of $A$:
	\begin{eqnarray*}
		A_{ii} &=&  \int_{i-1}^{i+1} (x+1) \phi'_i(x) \phi'_i(x) \\ 
		A_{ii} &=& \int_{{i-1}}^{i} (x+1) \phi_i'^2 dx + \int_{i}^{x_{i+1}} (x+1)\phi_i'^2 \\ 
		A_{ii} &=& \frac{1}{h_i^2}\int_{{i-1}}^{i} (x+1)  dx + \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx 
	\end{eqnarray*}

	\begin{eqnarray*}
		A_{i,i+1} &=& \int_{i}^{i+1} (x+1) \phi'_i(x) \phi'_{i+1}(x) \\ 
		A_{i,i+1} &=& \int_{i}^{i+1} (x+1) \phi_i' \phi_{i+1}' dx \\ %+ \int_{i}^{x_{i+1}} (x+1)\phi_i' \phi_{i+1}'  \\ 
		A_{i,i+1} &=& \frac{1}{h_i + h_{i+1}}\int_{i}^{i+1} (x+1)  dx %+ \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx 
	\end{eqnarray*}

	\begin{eqnarray*}
		A_{i,i-1} &=&  \int_{i-1}^{i} (x+1) \phi'_i(x) \phi'_{i-1}(x) \\ 
		A_{i,i-1} &=&  \int_{i-1}^{i} (x+1) \phi_i' \phi_{i-1}' dx \\%+ \int_{i}^{x_{i+1}} (x+1)\phi_i' \phi_{i-1}' \\ 
		A_{i,i-1} &=& \frac{1}{h_i + h_{i-1}}\int_{{i-1}}^{i} (x+1) dx % dx + \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx 
	\end{eqnarray*}
	
	These integrals are evaluated numerically in my program \textit{gfemddn} using Simpson's Rule.

	Solving for the forcing function $g$ in  \autoref{eq:scheme}
	\begin{eqnarray*}
		-\nabla \cdot \sigma(x) \nabla u(x) &=& g(x) \\
		-\frac{\partial}{\partial x} \cdot (x+1) \frac{\partial }{\partial x} \sin(x) &=& g(x) \\
		-\frac{\partial}{\partial x} \cdot (x+1) \cos(x) &=& g(x) \\
		 \sin(x) + x\sin(x) - \cos(x) &=& g(x) \\
		 \sin(x) (x+1) - \cos(x) &=& g(x) \\
	\end{eqnarray*}	 

	The load vector $g$ is then:
	\[
		g_i = <g, \phi_i> = \int_0^\pi \left(  \sin(x) (x+1) - \cos(x) \right) \phi_i
	\]
	In \textit{gfemddn} this integral is again evaluated using Simpson's Rule.

	\end{enumerate}