1 
% $Id$ 
2 

3 
\chapter{The First Steps} 
4 

5 
\section{Introduction} 
6 

7 
\subsection{Getting the software} 
8 

9 
\escript, \ESyS, all freely available. Where do people get \finley from? 
10 

11 
\begin{enumerate} 
12 
\item how to get the software 
13 
\item a few words about the general structure 
14 
\item installation 
15 
\end{enumerate} 
16 

17 
\section{How to solve a linear PDE} 
18 

19 
\begin{figure} 
20 
\centerline{\includegraphics[width=\figwidth]{FirstStepDomain}} 
21 
\caption{Domain $\Omega$ with outer normal field $n$.} 
22 
\label{fig:FirstSteps.1} 
23 
\end{figure} 
24 

25 
\begin{figure} 
26 
\centerline{\includegraphics[width=\figwidth]{FirstStepMesh}} 
27 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
28 
each element is a quadrilateral and described by four nodes, namely 
29 
the corner points. The solution is interpolated by a bilinear 
30 
polynomial.} 
31 
\label{fig:FirstSteps.2} 
32 
\end{figure} 
33 

34 
We want to solve the \index{partial differential equation}(\index{PDE}) 
35 
\begin{equation} 
36 
\Delta u + \alpha u =f 
37 
\label{eq:FirstSteps.1} 
38 
\end{equation} 
39 
for the solution $u$ on the domain $\Omega$. Here we assume that the 
40 
domain is the rectangle of length $1$ and height $2$, see 
41 
\fig{fig:FirstSteps.1}. $\Delta$ denotes the \index{Laplace 
42 
operator} which is defined by 
43 
\begin{equation} 
44 
\Delta u = (u\hackscore {,1})\hackscore{,1}+(u\hackscore{,2})\hackscore{,2} 
45 
\label{eq:FirstSteps.1.1} 
46 
\end{equation} 
47 
where for any function $w$ and any direction $i$ $u\hackscore{,i}$ 
48 
denotes the derivative of $w$ with respect to $i$. $\alpha$ is a 
49 
known constant (we will set $\alpha=10$) and $f$ is a given function 
50 
which may depend on the location in the domain. On the boundary of 
51 
the domain $\Omega$ the solution $u$ shall fulfill the socalled 
52 
homogeneous \index{Neumann boundary condition} 
53 
\begin{equation} 
54 
\frac{\partial u}{\partial n} = 0 
55 
\label{eq:FirstSteps.2} 
56 
\end{equation} 
57 
where $n=(n\hackscore1,n\hackscore2)$ denotes the outer normal field 
58 
of the domain, see Figure~\ref{fig:FirstSteps.1} and 
59 
\begin{equation} 
60 
\frac{\partial u}{\partial n} = n\hackscore1 u\hackscore{,1} + 
61 
n\hackscore2 u\hackscore{,2} 
62 
\label{eq:FirstSteps.2.1} 
63 
\end{equation} 
64 
denotes the normal derivative on the boundary. The partial 
65 
differential \eqn{eq:FirstSteps.1}) together with the 
66 
boundary condition~(\eqn{eq:FirstSteps.2}) forms a boundary value 
67 
problem (\index{BVP}) for unknown function $u$. 
68 

69 
In most cases, the BVP cannot be solved analytically and numerical 
70 
methods have to be used construct an approximation of the solution 
71 
$u$. Here we will use the \index{finite element method} 
72 
(\index{FEM}). The basic idea is to fill the domain with a set of 
73 
points, so called nodes. The solution is approximated by its values on 
74 
the nodes. Moreover, the domain is subdivide into small subdomain, 
75 
socalled elements. On each element the solution is represented by a 
76 
polynomial of a certain degree through its values at the nodes located 
77 
in the element. The nodes and its connection through elements is 
78 
called a \index{mesh}. \fig{fig:FirstSteps.2} shows an example 
79 
of a FEM mesh with four elements in the $x_0$ and for elements in the 
80 
$x_1$ direction over a rectangular domain. On more details we referring 
81 
to the literature, for instance \cite{X1,X2,X3}. 
82 

83 
\escript provides the class \linearPDE to define a 
84 
general linear, steady differential partial differential equation of 
85 
second order. We will discuss the most general form that can be 
86 
defined through \class{linearPDE} later. The components which are 
87 
relevant for us here is as follows: 
88 
\begin{equation} 
89 
\sum\hackscore{i,j=0}^k (A\hackscore{ij} 
90 
u\hackscore{,j})\hackscore{,i} + D u = Y 
91 
\label{eq:FirstSteps.3} 
92 
\end{equation} 
93 
In this form $D$ and $Y$ are scalars and $A$ is a $k \times k$ matrix 
94 
where $k$ denotes the spatial dimension (in our example we have 
95 
$k=2$). By comparing the template~(\ref{eq:FirstSteps.3}) with the 
96 
differential equation~(\ref{eq:FirstSteps.1}) we want to solve we can 
97 
immediately identify the appropriate values for $A$, $D$ and $Y$: 
98 
\begin{equation} 
99 
\begin{array}{lcccc} 
100 
A\hackscore{ij} & = & \delta\hackscore{ij}& =& 
101 
\left[ 
102 
\begin{array}{cc} 
103 
1 & 0\\ 
104 
0 & 1 
105 
\end{array} 
106 
\right] \\ 
107 
D & =& \alpha \\ 
108 
Y & = & f \\ 
109 
\end{array} 
110 
\label{eq:FirstSteps.3.1} 
111 
\end{equation} 
112 
When the PDE is defined via~(\eqn{eq:FirstSteps.3}), \class{linearPDE} 
113 
makes a particular assumptions about the boundary conditions: 
114 
\begin{equation} 
115 
\sum\hackscore{i,j=0}^k n\hackscore{i} A\hackscore{ij} 
116 
u\hackscore{,j} = 0 
117 
\label{eq:FirstSteps.4} 
118 
\end{equation} 
119 
Note that this boundary condition does not require extra information 
120 
as it only refers to the coefficient $A$ which already appears in the 
121 
PDE and so natural for it. Therefore this boundary condition is called 
122 
a \index{natural boundary condition}. 
123 

124 
We will set $\alpha=10$ and $f=10$ such that $u=1$ becomes the exact 
125 
solution of the boundary value problem. We make this very simple 
126 
choice to be able to test our program as we can compare the result 
127 
with the exact solution. Later we will set $\alpha$ and $f$ to 
128 
functions of their locations in the domain in which case we will not 
129 
be able to give an analytic solution. However, after testing our 
130 
program on this very simple case, we can be confident that it working 
131 
correctly before we apply it is a more complicated situation. 
132 

133 
This is the program to solve the boundary value problem: (Remember that 
134 
lines starting with '\#' are commend lines in Python) 
135 
%\verbatiminput{examples/FirstSteps1.py} 
136 
\begin{python} 
137 
# import ESyS and finley 
138 
from ESyS import * 
139 
import finley 
140 
# set a value for alpha: 
141 
alpha=10 
142 
# generate mesh: 
143 
mydomain=finley.Rectangle(n0=40,n1=20,l0=2.,l1=1.) 
144 
# generate a system: 
145 
mypde=linearPDE(A=[[1,0],[0,1]],D=alpha,Y=10,domain=mydomain) 
146 
# generate a test solution: 
147 
u=mypde.getSolution() 
148 
# calculate the error of the solution 
149 
error=u1. 
150 
print "norm of the approximation error is ",Lsup(error) 
151 
\end{python} 
152 
Line 2 import \escript and a few other tools for \ESyS. In Line 3 is 
153 
importing \finley which is used to solve the partial differential 
154 
equation. In line 7, a rectangular domain of length $l\hackscore 0=2$ 
155 
and height $l\hackscore 1=1$ is generated and subdivided in 
156 
$n\hackscore 0=40$ and $n\hackscore 1=20$ elements in $x\hackscore 0$ 
157 
and $x\hackscore 1$ direction, respectively. 
158 

159 
We are using a function of \finley. This determines later in the code 
160 
which solver for the PDE is actually being used to solve the PDE. 
161 

162 
The solution is done in three steps: 
163 
\begin{enumerate} 
164 
\item generate a finite element mesh subdividing the domain into elements. 
165 
\item assemble the system of linear equations $Mu=b$ from the BVP 
166 
\item solve the linear system to get $u$ 
167 
\end{enumerate} 
168 
The returned $u$ is given an approximation of the solution of the BVP 
169 
at the nodes of the finite element mesh. The quality of the 
170 
approximation depends on the size of the elements of the finite 
171 
element mesh: As smaller the element size the better the 
172 
approximation. In our example we know the solution of the BVP so we 
173 
can compare the returned approximation with the true solution. In 
174 
fact, as the true solution is simple, we can expect that the 
175 
approximation is exact. 
176 

177 
The first step imports the package \ESyS which includes among 
178 
others the module \finley. 
179 

180 
\section{A Time Dependent Problem} 
181 

182 
% \verbatiminput{exam/finley\hackscoretime.py} 
183 

184 
\section{With Dirichlet Conditions} 
185 

186 
% \verbatiminput{exam/finley\hackscoredirichlet.py} 
187 

188 
\section{Systems of PDEs} 
189 

190 
% \verbatiminput{exam/system\hackscoretime.py} 
191 

192 
\section{Explicit Schemes} 
193 
% \verbatiminput{exam/explicit\hackscoretime.py} 