Faculty of Science

Numerically solving the wave equation using the

finite element method

Bachelor’s Thesis

Jasper Everink
Mathematics
Supervisor:
Dr. Tristan van Leeuwen
Mathematical Institute

7th June 2018

CONTENTS i

Contents
0 Introduction 1
0.1 Wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
0.2 Finite element method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
0.3 This thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1 One-dimensional wave equation 3

1.1 Weak formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Nodes and Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Finite element solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 Time discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.6 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Irregular elements and higher-order approximation 11

2.1 Canonical element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Higher-order polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Matrix construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Numerical analysis and experiments 18

3.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Conservation of energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4 Two-dimensional wave equation 26

4.1 Higher-dimensional weak formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2 Elements and Nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.3 Canonical Elements and Basis function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.4 Finite Element Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5 Conclusion 36

A Implementation 37
A.1 One-dimensional wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
A.2 Irregular elements and higher-order approximation . . . . . . . . . . . . . . . . . . . . . . . . 37
A.3 Numerical analysis and experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
A.4 Two-dimensional wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

References I
0 INTRODUCTION 1

0 Introduction
0.1 Wave equation
The wave equation is a partial differential equation that is used in many field of physics. It is used to model
different kind of waves, for example sound waves and the oscillation of strings. The wave equation can be
defined in d-dimensions as
B2 u
c2 ∆u  0 with c ¡ 0,
B t2 (0.1.1)

where u is a function from Rd  R to R and ∆ is the Laplacian over the spatial coordinates defined as

∆u 
B2 u .
¸
d

n1
Bx2n (0.1.2)

A special case is the one-dimensional wave equation which can be more simply written as

B2 u
c2 B2 u  0 c ¡ 0.
B t2 B x2 with (0.1.3)

The function u can be the displacement of a string or membrane of a drum, the height of a water wave or
air pressure for sound waves, dependant on the field the equation is used in.

The constant c in the above equations defines the speed at which the wave ”moves”. To see this we can
look at the one dimensional wave equation for which, if we have no other conditions, upx, tq  f px ctq and
upx, tq  f px
ctq are solutions if the function f pη q is twice differentiable. These two solutions have a con-
stant shape which move in time with speed c to the left and right respectively. A general solution to the one-
dimensional wave equation was actually derived by d’Alembert and is given by upx, tq  f px
ctq g px ctq,
being the sum of a left and right moving wave.

In this thesis the focus will be on the one-and two-dimensional wave equations over a bounded domain D
where the spatial boundary of the domain B D is fixed at zero and we have an initial value in time. These
initial boundary conditions are specified as follows

upx, tq  0 for all x P BD and t P R,

upx, 0q  u0 pxq for all xPD and (0.1.4)
Bu px, 0q  u pxq x P D.
Bt 1 for all

Here, upx, tq represents the the displacement of a string (in one-dimension) or membrane of a drum (in
two-dimensions) that is held fixed at the boundary.

0.2 Finite element method

The finite element method is a numerical method generally used to solve differential equations with boundary
conditions. The concept of the finite element method is to divide the domain into finitely many smaller areas
called elements and approximate the solution of the differential equation on these elements using a suitable
set of basis functions. The result is a system of equations which can either be solved directly or by another
numerical method.

The finite element method has a few advantages over other numerical methods for solving differential equa-
tions (like the finite difference method). Firstly, the finite element method approximates the solution for
almost1 the whole domain, while most other methods only approximate the solution at a discrete set of
points. Secondly, with the finite element method we have more freedom over the discretization of the do-
main. We have the choice to increase precision at parts of the domain by using smaller elements and to
1 In chapter 4 we will see some cases where we cannot easily find an approximation for over the whole domain.
0 INTRODUCTION 2

decrease the precision at parts by using larger elements. Thirdly we can choose different basis functions on
the elements to increase or decrease the precision at different parts of the domain or choose basis functions
which better approximate the actual solution, for example by using trigonometric functions instead of poly-
nomials.

In this thesis we will use the finite element method for spatial discretization of the wave equation 0.1.1
and solve the resulting system of ordinary differential equations using a finite difference method. In the
one-dimensional case the elements will be intervals and we will use arbitrary degree polynomials for the
approximation. In the two-dimensional case the elements we will consider are triangular and we will use
linear approximation.

0.3 This thesis

Most of the practical side of this thesis is based on ”Finite Element Methods: A Practical Guide”[1], but
explained through the wave equation. The first two chapters of this thesis cover the one-dimensional case
and are based on chapter 2, 3 and 4 of that book[1]. In chapter 3 the numerical method for the one-
dimensional case is analysed and certain convergence and stability properties are discussed. Chapter 4
discusses a generalisation of the finite element method to the two-dimensional wave equation. The chapters
based on the book[1] use the basic steps as described in subsection 2.9 of that book, which can be summarized
as follows:
1 Derive a ”weaker” formulation of the differential equation.

2 Define the spatial discretization of the domain into elements and nodes.
3 Define the basis functions for the approximation.
4,5 Derive a system of equations.

6 Solve the system of equations.

The numerical methods derived in the chapters have also been implemented in Python and are used to make
most of the plots. An explanation of some important parts of the code can be found in appendix A and the
implementations can be found at www.github.com/jeverink/BachelorsThesis.
1 ONE-DIMENSIONAL WAVE EQUATION 3

1 One-dimensional wave equation

The first 4 sections of this chapter are based on chapters 2 and 3 of [1].

In this chapter a numerical method is derived for the wave equation in one dimension on the unit interval r0, 1s.

The one-dimensional wave equation is defined as follows:

B2 u
c2 B2 u  0 for c ¡ 0,
B t2 B x2 (1.0.1)

with initial boundary conditions

up0, tq  up1, tq  0 for all t P R,

upx, 0q  u0 pxq for all x P r0, 1s and (1.0.2)
Bu px, 0q  u pxq x P r0, 1s.
Bt 1 for all

1.1 Weak formulation

First we will derive a different but equivalent formulation for the one-dimensional wave equation defined in
1.0.1. Let v pxq be a differentiable functions such that v p0q  v p1q  0 and therefore zero on the boundaries
of the domain. By multiplying the wave equation with this function and integrating it over the unit interval
we get following expression »1 2
0
B u px, tq
c2 B2 u px, tq vpxqdx.
0 B t2 B x2 (1.1.1)

This expression can be rewritten, using integration by parts, to

»1   1 » 1
0
B 2
u
px, tqv pxqdx
c2 B u
p x, tqv pxq

Bu px, tq Bv px, tqdx
0 Bt Bx 0 Bx Bx
2
0
»1 2 »1
 BBtu2 px, tqvpxqdx c2 BBux px, tq BBxv px, tqdx
(1.1.2)
0 0
 apupx, tq, vpxqq,
»1 2 »1
with apupx, tq, v pxqq :
B u
px, tqvpxqdx c BBux px, tq BBxv px, tqdx.
2
0 B t 2
0
(1.1.3)

For the integrals above to exist we require certain properties of u, v and their derivatives.
For this we will use the L2 -norm on a domain D. For a function f this is defined by
»
||f ||2 : |f pxq|2 dx. (1.1.4)
D

Using this norm a function f is called square-integrable if the L2 -norm of f is finite or more specific

||f ||2 8. (1.1.5)

The set of all square-integrable functions is denoted by L2 . An important property of square-integrable

functions is that the integral of the product of two square-integrable functions over the domain D is also
finite. This property is important for the existence of the integrals in 1.1.2.

We want the functions u, v and their first derivatives with respect to x to be square-integrable. The set
of functions which are square-integrable and whose first derivative is square-integrable is called a first order
Sobelov space and can be denoted by H 1 . The set of functions which are also zero at the boundary of the
domain is denoted by H01 .
1 ONE-DIMENSIONAL WAVE EQUATION 4

We now want to find a function u which is in H01 for each moment in time t such that for every function v
in H01 the equations of 1.1.2 hold.
This gives rise to the weak formulation of our one-dimensional wave equation, which can now be stated as
follows

Definition 1.1.1 (Weak formulation). Find upx, tq such that for all t P R it holds that up, tq P H01 and for
all v P H01
apupx, tq, v pxqq  0, (1.1.6)
and the initial boundary conditions 1.0.2 holds.
This is called the weak formulation as the function u only has to hold in 1.1.6 with respect to certain functions
v called test functions.
Except for the derivation of the weak formulation above the rigorous connection between the weak formulation
and the ”strong formulation” 1.0.1 and 1.0.2 is beyond the scope of this thesis.

1.2 Nodes and Elements

We will now divide the interval r0, 1s into N elements of equal size separated by nodes.
The size of each element is given by h  N1 , therefore the nodes can be defined as

xi : pi
1qh for i  1, . . . , N 1, (1.2.1)

and the elements are defined as the intervals

Ei : rxi , xi 1 s for i  1, . . . , N. (1.2.2)

The nodes and elements are illustrated in the figure below.

Irregular elements are treated in Chapter 2.

Ei
1 Ei
x1 0 xi
1 xi xi 1 xN 1 1
h h

Figure 1.1: A regular partitioning of the unit interval r0, 1s in elements and nodes.

1.3 Basis functions

Instead of searching for a function u in H01 such that for all test functions v in H01 the equation of the weak
formulation hold, we are going to search in a finite-dimensional subspace of H 1 . This space will consist
of those functions in H 1 which are piecewise polynomial on each element. In this chapter we will restrict
ourselves to piecewise linear functions. We will denote this finite-dimensional set of functions as S 1 .

The set S 1 will be defined as the functions v that can be written as a linear combination of finitely many
basis functions and that are zero at the boundaries. We will take N 1 basis functions, one for each node,
and call these φi pxq for i  1, . . . , N 1. Then the functions v in S 1 can be written as
N¸1
v pxq  vi φi pxq, with vi PR and

i 1
(1.3.1)
v p0q  v p1q  0.

To get piecewise linear functions on the elements the basis functions are defined as follows.
For j  2, . . . , N we will define a basis function which are non-zero and linear on the elements Ej
1 and
1 ONE-DIMENSIONAL WAVE EQUATION 5

Ej and zero everywhere else. These basis functions can be written as

$ x
x
'
& h
j
1
xj
1 ¤ x ¤ xj
φj pxq  xj 1
x x ¤ x ¤ x , for j  2, . . . , N. (1.3.2)
j j 1
'
%
h
0 otherwise
For j  1, N 1 we define basis functions that lie on the boundary of the interval, such that they are non-zero
on the elements E1 and EN respectively and zero everywhere else. These functions can be written as
#
x2 x
x1 ¤ x ¤ x2
φ1 pxq  h ,
0 otherwise
# (1.3.3)
x xN
xN ¤ x ¤ xN
φN 1 pxq  0
h
otherwise
1
.

These basis functions are illustrated in figure 1.2.

An important property of these basis functions is that φi pxq is precisely one at node xi and zero on every
other node, therefore #
1 when i  j
φi p x j q 
0 when i  j
. (1.3.4)

This results in the property v pxi q  vi for i  1, . . . , N 1.

φ1 p x q φj p x q φN 1 pxq
1

0
x1 x2 xj
1 xj xj 1 xN xN 1
0 E1 E j
1 Ej EN 1

Figure 1.2: The piecewise linear basis functions φ1 , φj and φN 1 .

1.4 Finite element solution

Using the just defined basis functions we can define a finite element solution and finite element test functions
as follows
N¸1
U px, tq : Ui ptqφi pxq, (1.4.1)

i 1
and
N¸1
V pxq : Vi φi pxq. (1.4.2)

i 1

Instead of searching for a solution in H 1 we will search in S 1 . A finite element formulation can then be
formulated as follows
1 ONE-DIMENSIONAL WAVE EQUATION 6

Definition 1.4.1 (Finite element formulation). Find U px, tq such that for all t P R it holds that U p, tq P S 1
and for all V P S 1
apU px, tq, V pxqq  0, (1.4.3)
and the following initial boundary conditions hold

U p0, tq  U p1, tq  0 for all t P R,

U px, 0q  U0,0 pxq for all x P r0, 1s and (1.4.4)
BU px, 0q  U pxq x P r0, 1s.
Bt 0,1 for all

Note that because of the form U takes, it does not have to hold that U0 pxq equals u0 pxq and initially
we will choose U0 pxq and U1 pxq such that for each node xi it holds that Ui p0q  U0,0 pxi q  u0 pxi q and
BUi p0q  U0,1 pxi q  u1 pxi q. In chapter 3 we will use a more complex initial value to simplify the numerical
Bt
analysis.

For the finite element solution it is enough to test it only against the basis functions. This follows from the
linearity of apu, v q in the second variable and therefore it holds that
N¸1
0  apU px, tq, V pxqq  apU px, tq, Vi φi pxqq

i 1
(1.4.5)
N¸1
 Vi apU px, tq, φi pxqq.

i 1

Thus, if it holds for all of the basis functions φi , it also holds for all functions V pxq in S 1 .
For each basis function φi pxq that is not on the edge with i  2, . . . , N we can write

0  apU px, tq, φj pxqq

» 1 N¸1 » 1 N¸1
 B2 U ptqφ pxqφ pxq 2 B φ pxq B φ pxq
Ui ptq
0 
i 1
B t2 i i j dx c
0 i1
Bx i Bx j dx


N¸1
B2 U ptq » 1 φ pxqφ pxqdx 2
N¸1
Ui ptq
»1
B φ pxq B φ pxqdx

i 1
B t2 i 0 i j c
i1 0 B x
i
Bx j (1.4.6)
N¸1 N¸1
 : ptqT
U i i,j c2 Ui ptqSi,j ,

i 1 i 1 
»1 »1
with Ti,j : φi pxqφj pxqdx and Si,j :
B φ pxq B φ pxqdx.
0 0 Bx i Bx j
This results in N
1 linear ordinary differential equations. For the two basis functions left out we can use
the following additional equations
: p tq  0
N 1 ptq  0.
U and U : (1.4.7)
1

to make sure the boundary conditions hold.

By defining the vector U ptq : pU1 ptq, . . . , UN 1 ptqq we can use these N 1 linear differential equations to
write the following ordinary differential equation
: ptq
TU c2 SU ptq  0, (1.4.8)

with the coefficients of T and S being those derived in 1.4.6 except for the first and last row whose elements
are derived from 1.4.7.
The values of Ti,j and Si,j can easily be calculated by splitting their integrals over all the elements like
»1 N » xi
¸
f pxqdx  f pxqdx.
1

(1.4.9)
0 
i 1 xi
1 ONE-DIMENSIONAL WAVE EQUATION 7

As the basis functions only have a maximum of two elements on which they are non-zero, the values of
Ti,j and Si,j can be calculated by just integrating over the few elements on which both basis functions are
non-zero. This allows us to calculate the following coefficients
» xi » xi
 φi pxq2 dx φi pxq2 dx 
2
1

Ti,i h,
xi
1 xi 3
» xi
Ti,i
1  φi pxqφi
1 pxqdx 
1
6
h,
xi
1
» xi 
Bφi pxq 2 dx » x  Bφi pxq 2 (1.4.10)
 dx 
2
i 1

Si,i
xi
1 Bx x i
Bx h
and
» xi
Si,i
1  Bφi pxq Bφi
1 pxqdx 
1 .
xi
1 Bx Bx h

For all other values the basis functions do not overlap on non-zero elements and thus their product is zero.
A more general way of calculating these coefficients is given in chapter 2.

The result is a second order linear differential equation with constant coefficients 1.4.8. The solution to this
equation can than be used in the finite element solution to get a numerical solution to the one-dimensional
wave equation stated in the beginning of the chapter.

1.5 Time discretisation

In this section we derive two numerical methods to solve differential equation 1.4.8.

For the first method we will use a second order approximation for the second derivative given by

: ptq 
U pt dtq
2U ptq U pt
dtq
U , (1.5.1)
dt2
and substituting this in the differential equation results in

U pt dtq
2U ptq U pt
dtq
T c2 SU ptq  0. (1.5.2)
dt2

After some algebraic manipulation the following equation can be found

U pt dtq 
c2 dt2 T
1 SU ptq 2U ptq
U pt
dtq. (1.5.3)

As the initial values are U p0q and U9 p0q we can use a first order approximation to do the first step

U pdtq  U p0q dtU9 p0q, (1.5.4)

and use equation 1.5.3 for further iterations.

Notice that 1.5.3 requires the inverse of the sparse (or more specific tridiagonal) matrix T and its inverse is
generally not sparse. For a more efficient computation we can use

U pt dtq
2U ptq U pt
dtq
T
dt2

c2 SU ptq. (1.5.5)

By first calculating the right-hand side we get a linear system which we can solve for U pt dtq
2U ptq U pt
dtq
from which we can calculate U pt dtq. Using this method the inverse of a sparse matrix does not have to
dt2

be calculated, but instead a linear system has to be solved each iteration, which is in general more accurate
than a matrix inversion.
1 ONE-DIMENSIONAL WAVE EQUATION 8

A second method for solving equation 1.4.8 is to rewrite it as a first order linear differential equation. First
define W ptq : U9 ptq, then because W : ptq 
c2 T
1 SU ptq the differential equation is equivalent to the
9 ptq  U

initial value problem

# #
W9 ptq 
c2 T
1 SU ptq, U p0q  U0,0 and
U9 ptq  W ptq W p0q  U p0q  U0,1 .
with initial value (1.5.6)
9

Then by using a first order approximation

U pt dtq
U ptq
U9 ptq  , (1.5.7)
dt
we get the equations #
p
W t dt W t q
p q 
c2 T
1 SU ptq,
p q
p q
dt
U t dt U t
dt  W ptq, (1.5.8)

which can be rewritten to #

T p
U9 t dt U9 t q
p q 
c2 SU ptq,
dt
(1.5.9)
U pt dtq  U ptq dtU ptq.
9

For the first equation begin by solving the linear system then solve for U9 pt dtq and the second equation
can be calculated directly.

One advantage of the second method is that we do not have to use the initial value to the first step in time
and a second advantage will be clear in chapter 3 where we use the U ptq and U9 ptq to calculate the energy of
the solution.

1.6 Example
We will now use the second numerical method derived in the previous section to simulate a simple wave. As
an example let c  1 so the one-dimensional wave equation becomes

B2 u
B2 u  0.
B t 2 B x2 (1.6.1)

For this equation one of the solutions is given by

upx, tq  cosp2πtq sinp2πxq, (1.6.2)

which is periodic in time with period 1 and is called a standing wave.

The figures 1.3a and 1.3b show the simulation with time step 0.001 at different points in time for four and
respectively twenty elements. It can be seen in 1.3a that the numerical solution has a faster ”periodic”
movement than the actual solution. The word periodic is put between quotation marks as in 1.4 it can be
seen that the amplitude of the numerical solution increases over time.
In figure 1.3b it can be seen that the numerical method more accurately represents the actual solution, but
that it also moves faster than the actual solution. And just like with the four element example, the amplitude
increases over time as can be seen in 1.5
The behaviours just observed like the periodic movement, increase in amplitude and convergence will be
further analysed in chapter 3.
1 ONE-DIMENSIONAL WAVE EQUATION 9

(a) A standing wave simulated with 4 elements. (b) A standing wave simulated with 20 elements.

Figure 1.3: The standing wave upx, tq  cosp2πtq sinp2πxq simulated with a time step of 0.001 with 4 (left)
and 20 (right) elements.
1 ONE-DIMENSIONAL WAVE EQUATION 10

Figure 1.4: A standing wave simulated using 4 elements, showing that the amplitude of the numerical solutions
increased over time.

Figure 1.5: A standing wave simulated using 20 elements, showing that the amplitude of the numerical
solutions increased over time.
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 11

2 Irregular elements and higher-order approximation

The first 3 sections of this chapter are based on chapter 4 of [1].

In this chapter we will generalize the finite element method by first allowing irregular elements and then
higher order polynomial approximation.

We will divide an interval ra, bs with a b into N elements seperated by nodes. The nodes xi with i 
1, . . . , N 1 are chosen such that a  x1 x2 ... xN xN 1  b and the elements can than defined
the same as 1.2.2, namely
Ei : rxi , xi 1 s for i  1, . . . , N, (2.0.1)
and the size of element Ei is given by hi or more specific

hi : xi 1
xi . (2.0.2)

This is illustrated in figure 2.1;

Ei
1 Ei
x1 a xi
1 xi xi 1 xN 1 b
hi
1 hi

Figure 2.1: An irregular partitioning of the interval ra, bs in elements and nodes.

2.1 Canonical element

Instead of defining each basis function separately we will define a canonical element denoted Ec which is
defined as the unit interval r0, 1s. Then, for each element Ei there exists a transformation τi pxq such that
τi pEc q  Ei and the inverse of τi exists as illustrated in 2.2. The transformation τi is given by

τi pxq  xpxi 1
xi q xi  xhi xi . (2.1.1)

0 1
Ec

τi τj

Ei Ej
x1 xi xi 1 xj xj 1 xN 1
a b
hi hj

Figure 2.2: The transformations τi and τj from the canonical element Ec to the irregular elements Ei
and Ej .

On the canonical element we can define local basis functions

φl,1 pxq  1
x and
φl,2 pxq  x,
(2.1.2)
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 12

which are illustrated in figure 1.2.

1
φl,1 φl,2

0
Ec
0 1

Figure 2.3: The local linear basis functions φl,1 and φl,2 .

The global basis functions are then definable as

$
'
1
p p qq xi
1 ¤ x ¤ xi ,
&φl,2 τi
1 x
φi pxq  φl,1 pτi
1 pxqq xi ¤ x ¤ xi 1 , for i  1, . . . , N 1, (2.1.3)
'%
0 otherwise
with exceptions for basis functions corresponding to nodes on the boundary, like in chapter 1.

Using this definition the time and space matrix coefficients from 1.4.10 can be calculated using the following
change of variables for a functions f over an element Ek to the canonical element Ec
» xk
f pxqdx
1

xk
» τk 1 pq
 f pxqdx
τk 0 pq (2.1.4)
»1
 f pτk pxqqτ 1 pxqdxk
0
»1
hk f pτk pxqqdx.
0

If the function f pxq is of the the form g pτk
1 pxqq on element Ek (for example f pxq  φl,2 pτk
1 pxqqφl,1 pτk
1 pxqq
when calculating the coefficients Ti,j ) the above integral simplifies to
» xk »1
g pτk
1 pxqqdx  hk g pxqdx,
1

(2.1.5)
xk 0

which can be calculated easily when g pxq is a polynomial like is the case in this chapter. For example take
the regular grid (hi  h) of previous chapter, the coefficient Ti,i can be calculated like
» xi » xi
 φi pxq dx φi pxq2 dx
1
2
Ti,i
xi
1 xi
»1 »1
h φl,2 pxq2 dx h φl,1 pxq2 dx (2.1.6)
0 0
»1 »1
h x2 dx h p1
xq2 dx  h 32 .
0 0

or when using irregular elements

Ti,i  31 hi
1 1
3
hi . (2.1.7)
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 13

2.2 Higher-order polynomials

To increase the order of the polynomials we need more nodes. For dividing the interval ra, bs in N elements
and to approximate the solution on each element using a pth degree polynomial we need a total of N p 1
nodes. The nodes xi with i  1, . . . , N p 1 such that a  x1    xN 1  b can be split into two
groups: exterior and interior nodes.
The exterior nodes are xip 1 with i  0, . . . , N divide the interval into the elements such that

Ei : rxip p
1, x i 1 p 1q s and
for i  1, . . . , N. (2.2.1)
hi  xpi q
1 p 1
xip 1

For each element Ej there exist p
1 interior nodes xjp 1 i with i  1, . . . , p
1. For our purposes we will
place these interior nodes uniformly on the elements such that

: xjp  1, and i  1, . . . , p
1.
hj
xjp 1 i 1 i for j ... ,N (2.2.2)
p
The interior nodes can also be expressed on the canonical element The nodes and elements are illustrated in
figure 2.4.

Ei
1 Ei

x1 x3pi
1q 1 x3pi
1q 2 x3pi
1q 3 x3i 1 x3i 2 x3i 3 x3pi q

1 1 x3N 1
a b
hi
1 hi

Figure 2.4: An irregular grid for cubic approximation (p = 3).

We now want to define a basis function for each of the nodes. First we will define the local basis functions on
the canonical element. As there are p 1 exterior and interior nodes inside each element we need to define
p 1 local basis functions.
Let the local nodes xl,i be defined as

i
1
xl,i : for i  1, . . . , p 1. (2.2.3)
p
We want the corresponding local basis functions φl,i to be 1 on the local node xl,i and zero on all other nodes.
A suitable set of basis functions are Lagrange polynomials which result in the local basis functions
¹ x
xl,j
φl,i pxq :
¤¤
1 j p
x
1 l,i

xl,j . (2.2.4)

j i

For example, for quadratic approximation (p  2) there are three local basis functions

 px

21 qp
x 1
1q  2x2
3x
1
φl,1 1,
2
xpx
1q
φl,2  1 1 
4x
2

2 
2
4x and (2.2.5)

xpx
21 q
φl,3   2x2
x.
1 1 2

The quadratic and cubic cases are illustrated in 2.5 and 2.6.
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 14

φl,2 φl,3
1 1
φl,1 φl,2 φl,3 φl,1 φl,4

0 0
1 1 2
0 2 1 0 3 3 1

Figure 2.5: The quadratic local basis functions.

Figure 2.6: The cubic local basis functions.

For exterior nodes we can now define the corresponding global basis functions as
$
'
&φl,p pτi

11 pxqq
1 xi
1 ¤ x ¤ xi ,
φip 1 pxq  φl,1 pτi
1 pxqq xi ¤ x ¤ xi 1 , i  2, . . . , N, (2.2.6)
'%
0 otherwise

with the same boundary exceptions as earlier.

For interior nodes the global functions can be defined as
#
φl,i pτj
1 pxqq xj ¤ x ¤ xj
φjp 1 i pxq  0 otherwise
1,
j  2, ... ,N and i  2, . . . , p. (2.2.7)

Notice that the global basis functions only take non-zero value on the elements that the corresponding node
touches and again that for all i  1, . . . , N p 1 it holds that φi pxi q  1.

2.3 Matrix construction

³b B B ³b
A simple way of calculating the coefficients Ti,j
³b
 a
φi pxqφj pxqdx and Si,j 
Bx°φi pxq³Bx φj pxqdx is by a
N
rewriting the integrals a as a summation of the integrals over all the elements like i1 Ei and then use
the definition of the global basis functions and a change of variables to the canonical elements to simplify
the calculation of each integral. We can reduce the amount of integrals by not evaluating those for which we
know that the product of the functions is zero.

We can also calculate the coefficient the other way around. We know that the product of two basis functions
(and the product of their first derivatives to x) is non-zero on an element if and only if both corresponding
nodes touch the element (either as an external or as an internal node), therefore there are pp 1q2 pairs of
basis functions on each element whose product is non-zero on that element. These pairs are more specifically
pairs of transformed local basis functions on that element. Using this we can start with all values Ti,j and
Si,j zero and then iterate over all the element and add the integral of the non-zero pairs of that element to
the corresponding coefficients. Both methods (as example for Ti,j ) are summarized bellow:

Per coefficient calculation

All T[ i , j ] are 0 ;
For each c o e f f i c i e n t T [ i , j ] :


I f b a s i s f u n c t i o n s i and j have non z e r o o v e r l a p :
For each e l e m e n t E [ n ] :
T [ i , j ] += I n t e g r a t e b a s i s f u n c t i o n s i and j o v e r e l e m e n t E [ n ] ;
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 15

Per element calculation

Let p be t h e d e g r e e o f p o l y n o m i a l ;
All T[ i , j ] are 0 ;
For each e l e m e n t E [ n ] :
For each l o c a l b a s i s f u n c t i o n i :
For each l o c a l b a s i s f u n c t i o n j :
T [ ( p∗n)+ i , ( p∗m)+ j ] += I n t e g r a t e t r a n s f o r m e d l o c a l b a s i s f u n c t i o n s i and j o v e r
element E[ n ] ;

2.4 Examples
For the following examples we will choose the constant c  1, the domain will be the interval r0, 2s and the
time step will again be 0.001. The interval r0, 2s will be divided into elements of three different sizes, where
the smaller elements are on the left size and the largest element is on the right. This division is illustrated
by the green dots in 2.7 and 2.8.

For the first example we will use the same solution as in chapter 1 1.6.2 which is given by

upx, tq  cosp2πtq sinp2πxq. (2.4.1)

In 2.7a we can see that for first degree approximation we can see that the finer elements can better approx-
imate larger wave, while this is more difficult for the coarser elements. It is even possible that a part of
the wave is completely removed when the nodes of an element lie on the zeros of the initial value. As time
increases the wave becomes inaccurate quickly with respect to the actual solution.

If we increase the degree of approximation to four we get 2.7b. The increase in degree implies an increase
in the amount of nodes and smoothness and thus it initial approximation can better approximate the actual
initial value. It is important to notice when comparing 2.7a and 2.7b that even though the solution starts
better approximating the actual solution when the degree of polynomials is higher the finer approximation
gets unstable faster. This might be explained by the increment of calculations with relatively small values,
which can easily introduce rounding errors.

As a second example we will use the sum of two standing waves or more specifically

upx, tq  cosp2πtq sinp2πxq cosp4πtq sinp4πxq. (2.4.2)

This example is plotted for different degrees of approximation in 2.8. In this example we see better that in-
creasing the degree of the polynomials improves the numerical solution, but again, making the approximations
too fine can introduce instability as can be seen in 2.8d.
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 16

(a) A first degree polynomial approximation. (b) A fourth degree polynomial approximation.

Figure 2.7: The wave upx, tq  cosp2πtq sinp2πxq simulated using an irregular grid with first (left) and fourth
(right) degree polynomials.
2 IRREGULAR ELEMENTS AND HIGHER-ORDER APPROXIMATION 17

(a) A first degree polynomial approximation. (b) A second degree polynomial approximation.

(c) A third degree polynomial approximation. (d) A fourth degree polynomial approximation.

Figure 2.8: The wave upx, tq  cosp2πtq sinp2πxq cosp4πtq sinp4πxq simulated using an irregular grid with
first, second, third and fourth degree polynomials.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 18

3 Numerical analysis and experiments

In this chapter we will analyse some phenomena we have seen in the examples of the previous chapters. First
we will look at the convergence of the numerical solution, then analyse the conservation of energy of solutions
and finally we will take a look at the dispersion of waves.

3.1 Convergence
Let u be the solution of the one-dimensional wave equation 1.0.1, U be the solution of the finite element for-
mulation 1.4.1 and 1.4.8 (after spatial discretization) and let Û be the final numerical solution (after temporal
discretization). We will work on the unit interval r0, 1s with arbitrary elements of size hi and approximating
the solution using pth degree polynomials

To analyse the convergence we want to find an upper bound for the total error u
Û with respect to a certain
norm. We can split this problem in two parts by using the triangle inequality resulting in
   

u
Û  ¤ ||u
U || 
U
Û  , (3.1.1)

therefore bounding the total error by the sum of the spatial and temporal discretization error.

First we will look at the spatial discretization error. To simplify the analysis we will require the following
extra constraints for all test functions V
»1
BV BU0,0 dx  » 1 BV Bu0 dx
0 Bx Bx 0 Bx Bx
and
»1 »1 (3.1.2)
BV B U0,1 dx  BV Bu1 dx.
2

0 B x B xB t Bx Bx 0

Just like for the finite element formulation it is enough to test the constraints only for the basis functions,
which implies they also hold for all test functions. These constraints are to make sure the initial value of
the discretization is close enough to the true initial value. For the finite element formulation discussed in
previous chapters the above constraints become two linear systems which can be solved to get initial values
for the finite element solution.

Before defining the norm we will use to find an upper bound of the error, we will first define the spatial part
of the norm as follows d
»1 
pf pxqq2 BBfx pxq dx,
2
||f pq||H1 : (3.1.3)
0

which can be seen as a form of energy. Then, we can combine this norm with a supremum norm to get the
following norm for a set R „ R¥0 :

||f ||L8 pR,H q : sup ||f p, tq||H . (3.1.4)

1
P
t R
1

A useful property we will use is that if R1 „ R, then ||f ||L8 pR1 ,H1 q ¤ ||f ||L8 pR,H q , because we can never get
1
a bigger value in a supremum when removing elements.

Because we discretize time we only want to use the norm on the times we actually calculate and define the
set R of all steps in time starting at 0, with steps of time dt up until a constant T . We will use the norm
above with this set R to find an upper bound for the total error.

For the spatial discretization error we will use the following theorem.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 19

Theorem 3.1.1 (Theorem 2.1 from [2]). Let u be the solution of the one-dimensional wave equation, U be
the solution of the finite element formulation such that the constraints 3.1.2 hold, h be the size of the biggest
element and T ¡ 0 an arbitrary constant, then there exists a constant C ¡ 0 such that

||u
U ||L8 pr0,T s,H q 1 Chp . (3.1.5)

Because R € r0, T s, it follows directly that

||u
U ||L8 pR,H q ¤ ||u
U ||L8 pr0,T s,H q

1 1 Chp , (3.1.6)

with the constant C from the theorem above. This implies that the spatial discretization up until time T
converges to zero when the element sizes go to zero and it converges faster with higher order polynomials.
This can be written as ||u
U ||L8 pR,H 1 q  Ophp q.

Now we will look at the temporal discretization error. In chapter one we derived a method for solving the
:̄
equation 1.4.8 T U ptq c2 S Ū  0 using the equivalent formulation
#
:̄
U ptq
c2 T
1 S Ū ptq, , (3.1.7)
U ptq  U ptq
9̄ 9̄

°n
where Ū ptq is a vector of functions such that Û px, tq  i1 Ūi ptqφi pxq. We will focus on this second method
instead of the first method we derived, because this one will also be important in the next section.

The second method 1.5.9 is given by

$

&T U9̄ pt q
p q
dt U9̄ t

c2 S Ū ptq,
dt
(3.1.8)
%Ū t p dtq  Ū ptq dtU ptq, 9̄

or can be written as #
U9̄ pt dtq 
c2 dtT
1 S Ū ptq U9̄ ptq,
(3.1.9)
Ū pt dtq  Ū ptq dtU ptq. 9̄

This method was derived by using a truncated multidimensional version of the following forward difference
equation[3]
f pt dtq
f ptq dt d2 f
df
dt
ptq 
dt

2 dt2 pξptqq, ξptq P rt, t dts. (3.1.10)

If f is a vector of functions (f ptq : pf1 ptq, . . . , fn ptqq) and we apply the above formula to each function of
the vector, we get the following vector valued formula

f pt dtq
f ptq dt
df
dt
p tq 
dt

2 F ptq, with
" 2 * (3.1.11)
F ptq  pξi ptqq and ξi ptq P rt, t dts.
d fi
dt2 i1,...,n

Placing this formula in 3.1.7 for U9̄ we get

#
p
U9 t dt U9 t q
p q
dt F ptq

c2 T
1 SU ptq, , or
p q
p q
dt Gptq
dt 2
U t dt U t
dt 2  U ptq, 9

# (3.1.12)
U9 pt dtq 
c2 dtT
1 SU ptq U ptq dt2 F ptq,
2
9

U pt dtq  dtU ptq U ptq dt2 Gptq,

. 2
9
3 NUMERICAL ANALYSIS AND EXPERIMENTS 20

This results in the following errors

# #
U9 pt dtq
U9̄ pt dtq

dt2 F ptq, or U px, t dtq
U px, t dtq 
dt2 F px, tq,
2 9̂ 2
9
(3.1.13)
U pt dtq
Ū pt 
dt2 Gptq,
dtq U px, t dtq
Û px, t dtq 
dt2 Gpx, tq.
2 2

By doing this for n steps such that dt  Tn we get

#
U px, t ndtq
U px, t ndtq  n dt2 C1  T dt C1 ,
9 9̂
with C1 and C2 constants. (3.1.14)
U px, t dtq
Û px, t dtq  n dt2 C2  T dt C2 ,

This results in the global error  


U
Û L8 pR,H q  OpT dtq. (3.1.15)
1

The result is that, if we let T be a constant as has been done in Theorem 3.1.1, we get the total error
 

u
Û L8 pR,H q  Ophp dtq, (3.1.16)
1

which implies that, on the interval r0, T s, the error converges to zero if both h and dt tend to zero and
the convergence speed of the spatial discretization error can be improved by increasing the degree p of
polynomials.

3.2 Conservation of energy

In the previous section we have seen a few parameters that are related to the convergence of the numerical
solution, namely the size h of the biggest element, the time step dt and the degree p of the polynomials In
this section we will use the concept of energy conservation to see how these parameter have an impact on the
convergence and stability of the numerical solution.

The concept of conservation of energy is used all over physics and states most of the time that the total
amount of energy in a closed systems is constant. This conservation can also hold for the one-dimensional
wave equation and can be formulated as
Theorem 3.2.1 (Theorem and prove based on theorem 5.4 in [4]). For a solution u of the one-dimensional
wave equation there exists a constant C ¡ 0, independent of time t, such that

B p q2
 u


 u
c 
2 B p q2  C,
B
 t , t  x
,t
B  (3.2.1)

where |||| is the L2 -norm over the domain D.

Proof. Let v pxq  BBut px, tq be the test function in the weak formulation of the one-dimensional wave equation.
This gives
»
0
B2 upx, tq Bupx, tq dx c2 » Bupx, tq B Bupx, tq
D Bt2 Bt D Bx Bx Bt
» »
 12 BBt BupBx,t tq dx c2 21 BBt BuBpx,x tq dx
2 2

(3.2.2)
D D
 2  2
1 B  B u  2  B u
 
 2 Bt  Bt p, tq c  Bx p, tq ,


which implies that 

 u
 B p q2 
 u
c 
2 B p q2  C P R.
B
 t , t  B
x
,t  with C (3.2.3)
3 NUMERICAL ANALYSIS AND EXPERIMENTS 21

Notice that the proof of theorem 3.2.1 uses the weak formulation, thus when using v pxq  U9 px, tq in the finite
element formulation 1.4.3 it follows that

B p q2
 U

 U
c 
2 B p q2  C P R.
B
 t , t  B
x
,t  with C (3.2.4)

This implies conservation of energy for the finite element solution of the wave equation.
 2
From the above formulas we can define an energy function as E pu, tq : ||u9 p, tq|| c2  BBux p, tq , such that
2

for a solution u of the weak or finite element formulation there exists a constant C ¡ 0 such that E pu, tq  C
for all t.

The first error we can notice is that the initial value of the wave equation upx, 0q and BBut px, 0q and its discret-
ized value for the finite element formulation U px, 0q and BBUt px, 0q do not necessarily have the same energy,
or more specific, E pu, 0q and E pU, 0q do not have to be equal. But if we make the elements smaller and use
higher order polynomials for approximation, we can better approximate the original initial value with U and
BU and therefore better approximate the energy.
Bt
To analyse what will happen to the energy when we discretize the system of ordinary differential equations
1.4.8 we can use a °
more efficient way of calculating the energy. We can get this by putting the finite element
solution U px, tq  i1 Ui ptqφi pxq in the energy function, resulting in the following equations:
n



E pU, tq  U9 p, tq

2
  U
c2 
B p, tq2
Bx 
» ¸ » ¸
B
n 2 n 2

 Ui ptqφi pxq
9 dx c2
Ui ptq φi pxq
B dx
D 
i 1 D i 1
x
» » (3.2.5)

¸
n ¸ n
Ui ptqUj ptqφi pxqφj pxqdx 2
¸
n ¸
n
B B
Ui ptqUj ptq φi pxq φj pxqdx,
B B
9 9 c thus
 
D i 1j 1 D i1 j 1 x x
¸ ¸
n n ¸
n ¸
n
E pU, tq  U9 i ptqU9 j ptqTi,j c2 Ui ptqUj ptqSi,j
 
i 1j 1  
i 1j 1

where Ti,j and Si,j are defined as in 1.4.6. This is one of the main reasons why we introduced two numerical
methods for solving the ODE in chapter 1. The first method was more easily derived, but does not give us
concrete values for U9 i ptq, while the second method does give us these values. Combining that with the fact
that we already had to calculate most (if not all) values Ti,j and Si,j to construct the matrices of the ODE,
means that we have all the values to calculate the energy.

Because we know that real solutions of the wave equation and finite element formulation need to have a
constant energy, we can use the energy function as a measure of error over time when we have all the values
Ui and U9 i at a moment in time, like for the second method. We can even use this method this analyse
numerical solutions for which we don’t have the real solutions.

The first place where we saw something that we can analyse through energy is the increase of amplitude
over time we could see in figures 1.4 and 1.5 in chapter 1. If we take the same example solution and use ten
regularly spaces elements, we can see the energy using different time steps in figures 3.1a and 3.1b.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 22

(a) Energy of the example over a short interval for differ- (b) Energy of the example over a long interval for different
ent time steps. time steps.

Figure 3.1: Energy of the example for different time steps.

In the figures we can see that the energy increases over time, but that the amount of increase decreases with
the time step. This increase in energy can be seen through the slowly increasing amplitude of the numerical
solution. While in the beginning this energy increase is relatively unstable (the energy increases, but not
exponentially), we can see in the logarithmic plot 3.1b that after some time there is a huge change in beha-
viour. After some time the energy starts to increase exponentially, but this exploding change happens later
and is less steep when the time step decreases. We can actually see the beginning of this explosive behaviour
on the bottom right of figure 1.5, where the solution becomes slightly less sinusoidal and the end result can
be seen in 2.7a and 2.8c where the numerical solution suddenly turns chaotic.

This is also why we used a norm that only uses time up until a finite value T instead of all positive numbers.
For finite T the energy is bounded, but if we remove this restriction, it might happen that the energy diverges
to infinite and the solution will not converge.

Now let us take a look at a more complex situation where we do not know the real solution. For this we
will use an example where the initial value consists of two Gaussian peaks, like in figure 3.2, and where the
derivative in time is initially zero everywhere.

Figure 3.2: Two Gaussian peaks as the initial value for the example.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 23

We can see the energy over time for different time steps for 10 and 20 regular elements in figure 3.3 and 3.4
respectively. In these figures the instability becomes more clear and we can even notice that smaller element
like for figure 3.4 is more unstable than in 3.3. Because in both cases the energy increases less rapidly with
smaller times steps we need both the time step and the element size to decrease sufficiently fast to not increase
the instability.

(a) Energy of the example over a short interval for differ- (b) Energy of the example over a long interval for different
ent time steps. time steps.

Figure 3.3: The energy of the example simulated using 10 regular elements.

(a) Energy of the example over a short interval for differ- (b) Energy of the example over a long interval for different
ent time steps. time steps.

Figure 3.4: The energy of the example simulated using 20 regular elements.

Thus, we can conclude from the analysis of the energy that we require both the decrease in size of elements
and of the time step for the convergence of the numerical solution. The size of the elements and the degree
of the polynomials for approximation mainly contributes to the convergence of the energy to that of the real
energy and the decrease of the time step mainly contributes to decreasing the deviation from the constant
energy, but we need to decrease both values sufficiently fast to not increase the instability.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 24

3.3 Dispersion
One phenomenon we have not touched on yet is the different periodic movement a numerical solution might
make compared to the actual solution, like could be seen in figure 1.3. We will analyse this by searching a
relation between the spatial and temporal frequencies of waves.

Not every wave is a solution of the wave equation. To see this we will take a look at solutions of the form[5]

upx, tq  eipωt kx q  cospωt kxq i sinpωt kxq, (3.3.1)

where ω is the (temporal) frequency and k is called the wave number or spatial frequency. When we put this
solution in the wave equation, we can derive a relation as follows:

B2 eipωt kx q
c2 B2 eipωt kx q  0,
Bt2 Bx2

ω 2 c2 k 2 eipωt kx q  0,
(3.3.2)
ω 2
c k ,
2 2

ω  ck.
This relation between the frequency and the wave number is called the dispersion relation of the wave equa-
tion. We can already see in figure 1.3 that the same relation does not hold for the numerical method. The
problem we are going to analyse now is if there is such a relation for our numerical method and if there is,
how does it relate to the dispersion relation derived above.

The first step is to discretize the°solution 3.3.1. If we have n nodes and basis functions then the finite element
solution is written as U px, tq  i1 Ui ptqφi pxq, where Ui ptq is the value at the node xi at time t. From this
n

we get the discretization

Ui ptq  eipωt kxi q . (3.3.3)
The numerical method 1.5.9 we used could be written as
#

T p q
p q
U9 t dt U9 t

c2 SU ptq,
dt
(3.3.4)
U pt dtq  U ptq dtU ptq, 9

and if we put the just derived discretization in there we get the following equation
#

T iωeiωdt U t p q
iωU ptq 
c2 SU ptq,
dt
(3.3.5)
e iωdt
U ptq  p1 dtiωq U ptq.
Let us first take a look at the second equation in 3.3.5. This equation has to hold for all t and because
Ui ptq  eipωt kxi q  0 has no solution it implies that eiωdt  1 idtω. If we rewrite it using Euler’s formula
we get cospωdtq i sinpωdtq  1 iωdt. From this equation we can conclude that ω is either 0 or complex
with non-zero imaginary part, because if ω would be real it must hold that sinpωdtq  ωdt, which can only
hold for real numbers if ωdt  0 and if we let dt be zero we have no use for the numerical method. If we
choose ω to be zero then there is no temporal movement in the solution, so the only non-trivial option is for
ω to be complex with Impω q  0.

But what happens when we write ω as Repω q iImpω q with Impω q  0. If we put this in the discretized
solution we can rewrite it as follows

Ui ptq  eipRepωqt p q
iIm ω t kxi q  e
Impωqt eipRepωqt kxi q, (3.3.6)

then this equation has two possible scenarios (already excluding Impω q  0), either Impω q ¡ 0 in which case
the solution will converge to zero or Impω q 0 in which case the solution will diverge to infinite. Therefore,
if there exists such a numerical solution, it will always either diverge to infinity or dim out to zero.
3 NUMERICAL ANALYSIS AND EXPERIMENTS 25

Let us now take a look at the first equation in 3.3.5. We can rewrite this equation as

iω eiωdt
1
T U ptq 
c2 SU ptq. (3.3.7)
dt
As an example, we will use the situation as in chapter 1, where the domain is the unit interval r0, 1s, which we
divide in N equally sized elements of size h  N1 such that for j  1, . . . , N 1 there is a node xj  pj
1qh.
In this case we already calculated the values for T and S. If we combine this with equation 3.3.7 we get for
j  1, . . . , N the following non-trivial equations for which we can derive a dispersion relation as follows
 
iω eiωdt
1
hUj
1 ptq hUj ptq 1 ptq 
c
h1 Uj
1 ptq Uj ptq
Uj ptq
1 2 1 2 2 1
hUj 1 ,
dt 6 3 6 h h
 
iω eiωdt
1 1
ikh
dt 6
he
2
3
h
1 ikh ipωt
6
he e kxj q 
c2
h1 e
ikh 2
h

h1 eihk eipωt kxj q,
  (3.3.8)
iω eiωdt
1 2   
e
ikh 
c
e
ikh
1 ikh 2 2 ihk
h e e 2 ,
dt 6 3

iω eiωdt
1 h2 pcospkhq 2q  c2 dt pcospkhq
1q .
1
6
This results in the dispersion relation
 6c2 cospkhq
1
iω iωdt
dt
e
1  2
h cospkhq 2
. (3.3.9)

This relation approximates the dispersion relation of the wave equation ω  ck for small values of ω and
k. To notice this we will approximate equation 3.3.9. We will use the the first few terms of the Taylor
expansions of ex and cos pxq
1
cospxq 2 , which are given as

ex 1 x Opx2 q, (3.3.10)

and
cospxq
1 2

cospxq 2

x6 Opx3 q. (3.3.11)

Filling these approximations in the dispersion relation results in

ω2  c2 k2 , (3.3.12)

for small values of ωdt and kh.

Therefore we might expect the numerical solution to resemble the real solution when ωdt and kh are small.
This gives another reason why we need decrease both dt and h sufficiently to improve the accuracy of the
dispersion relation.

This chapter we have seen the importance of element sizes and the time step with regard to the accuracy of
the solution and the requirement to decrease values sufficiently for both convergence and stability.
4 TWO-DIMENSIONAL WAVE EQUATION 26

4 Two-dimensional wave equation

In this chapter we will numerically solve the two-dimensional wave equation using the same techniques of the
first two chapters. The two-dimensional wave equation over a domain D € R2 is defined as

B2 u
c2 ∆u  0 for c ¡ 0,
Bt2 (4.0.1)

with initial boundary conditions

upx, y, tq  0 for all px, yq P BD,

upx, y, 0q  u0 px, y q for all px, yq P D and (4.0.2)
∇upx, y, 0q  u1 px, y q for all px, yq P D.
The derivation is almost the same as in previous chapters, but the introduction of the second spatial dimension
creates a few extra difficulties, but luckily the resulting ordinary differential equation is of the same form as
in the one-dimensional case. In this chapter we will mainly focus on the differences and extra difficulties of
solving the two-dimensional wave equations compared to the one-dimensional wave equation.

4.1 Higher-dimensional weak formulation

First we will need a weak formulation of the two-dimensional wave equation. Because it is easier, we will
actually derive the weak formulation for the wave equation in arbitrary dimensions. For this we will need a
generalization of integration by parts called Green’s first identity[6], which is formulated as
» » »
ψ∆φdV ∇ψ  ∇φdV  ψ p∇φ  nqdV, (4.1.1)
D D BD
where n is the normal of the hypersurface B D.

Just like for the one-dimensional case we will multiply the wave equation with a test function v px, y q, which is
zero at the boundary of the domain, and integrate over the whole domain. This results in the follow equation:
»
0
B2 u px, y, tqvpx, yqdA
c2 » ∆upx, y, tqvpx, yqdA
D Bt
2
D
» » »
B 2
 Bt2 px, y, tqvpx, yqdA
c
u 2
v px, y q p∇upx, y, tq  nq dA
∇upx, y, tq  ∇v px, y qdA
»
D
»
B D D

 BBtu2 px, y, tqvpx, yqdA c2 ∇upx, y, tq  ∇vpx, yqdA

2
(4.1.2)
D D
 apupx, y, tq, vpx, yqq,
»
with apupx, y, tq, v px, y qq :
B2 u px, y, tqvpx, yqdA c2 » ∇upx, y, tq  ∇vpx, yqdA.
D B t2 D

We will not discuss the the requirements on upx, y, tq and v px, y q in depth, except for the requirement that
they both need to be square integrable on the domain such that the integrals above exist.

The result is the two-dimensional weak formulation, where we search for a function upx, y, tq, such that for
all test functions v px, y q it holds that apupx, y, tq, v px, y qq  0 and that the conditions 4.0.2 hold.

The space which contains u and the test functions is infinite-dimensional, therefore we will search in a finite-
dimensional subspace spanned by global basis functions φi px, y q. The result is the finite element formulation
where we search for the finite element solution
¸
n
U px, y, tq : Ui ptqφi px, y q, (4.1.3)

i 1
4 TWO-DIMENSIONAL WAVE EQUATION 27

such that for all finite element test functions

¸
n
V px, y q : Vi φi px, y q, (4.1.4)

i 1

it holds that
apU px, y, tq, V px, y qq  0, (4.1.5)
and because of the bilinearity of ap, q it is enough to only test U px, y, tq against the global basis functions
φi px, y q.

Before we can define the actual global basis functions we first need to define the elements we are going to
approximate the solution on.

4.2 Elements and Nodes

Let there be N elements and n nodes. The nodes are xi P R2 for i  1, . . . , n and the elements Ej € R2 for
j  1, . . . , N are defined by three nodes for which we will use the indices Ej,i such that the three corners
of element Ej are xEj,1 , xEj,2 and xEj,3 , which we will from now on write as xj,1 , xj,2 and xj,3 respectively.
It is important that all elements are orientated the same way, thus all corners are either indexed clockwise
or counter-clockwise. Not orientating the elements the same can lead to unexpected results. The figures
4.1a and 4.1b are examples of discretization of a square domain and partial discretization of a quarter circle
domain respectively.

x1 x2 x1
x2

E1
E1
E2
E2
x3

x3 x4 x4
(a) Discretization of a square (b) Partial discretization of a quarter circle

These examples show one of the problems that arises when we go from one to two dimensions. In one di-
mension, every connected domain is an interval which we can always fully subdivide in more intervals, but in
two dimensions we cannot always fully subdivide a connected domain in finitely many triangles like in figure
4.1b. We can however use more elements to better approximate the domain or use special curved elements,
but for simplicity we will only discuss triangular elements in this thesis.

Another problem we get when going to a higher dimension is that subdividing elements does not necessarily
increase accuracy. For example, in the case of a triangular domain we can continuously subdivide the elements
like in figure 4.2. The result is many stretched out elements and the approximation between the extremities of
the elements does not improve. Therefore, while in one dimension an important measure for the convergence
was the length of the intervals, in two (and higher) dimensions this generalizes to the diameter (furthest
distance between points) of the elements.
4 TWO-DIMENSIONAL WAVE EQUATION 28

(a) One triangle (b) Two triangles (c) four triangles (d) eight triangles

Figure 4.2: Repeated subdivision

4.3 Canonical Elements and Basis function

Just like in chapter 2 we will make sure that each element is just a transformation from a simpler canonical
element. As all our elements are triangles, we will define our element Ec as the triangle spanned by the points
p0, 0q, p1, 0q and p0, 1q as illustrated in figure 4.3.
p0, 1q
xl.3

xl,1 xl,2
p0, 0q p1, 0q
Figure 4.3: Triangular canonical element Ec

We can then get each triangular element spanned by the points xj,1 , xj,2 and xj,3 using the transformation
τj given by
τi pX, Y q  p1
X
Y qxj,1 Xxj,2 Y xj,3 , (4.3.1)
which has an inverse if the triangle is not degenerate (the three points do not lie on the same line). The
transformation is illustrated in figure 4.4.

xl,3 x2
x3

Ei
τi

Ec
x1

xl,1 xl,2

Figure 4.4: Transformation τi from the canonical element Ec to the element Ei .

Before defining the global basis functions we have to first define the local basis functions on the canonical
element. We again want the local basis functions to be linear on the canonical element, one on their respective
4 TWO-DIMENSIONAL WAVE EQUATION 29

point and zero on all other nodes. These basis functions can simply be defined as

φl,1 px, y q  1
x
y,
φl,2 px, y q  x and (4.3.2)
φl,3 px, y q  y,

and are illustrated in figure 4.5.

φl,1 px, y q φl,2 px, y q φl,3 px, y q

1 1 1

0 0 0
1 1 1
y y y

x 1 x 1 x 1
(a) Local basis function φl,1 px, y q (b) Local basis function φl,2 px, y q (c) Local basis function φl,3 px, y q

Figure 4.5: two-dimensional local basis functions

We now have everything to define the global basis functions. We will construct a single global basis function
for each node xi as follows. If the node xi is a corner point of an element Ej such that xi  xj,k , therefore
xi is the k th corner point of Ej , then the corresponding basis function φi px, y q is given on that element by

φi px, y q  φl,k pτj
1 px, y qq, (4.3.3)

and is zero on all elements the node is not a corner point of or is not covered by an element.

The definition of the basis functions does give rise to a problem called hanging nodes. In figure 4.6a we
see node xh , which is at the corner of the elements E1 and E2 , but lies on an edge of E3 , therefore the
corresponding basis function φh is one at xh for the elements E1 and E2 , but is zero on the whole of E3 ,
causing a discontinuity on the edges connecting E1 and E2 to E3 . If the value Uh ptq is non-zero, then this
discontinuity will also be in the finite element solution. Because we do not want this discontinuity, we also
do not want these hanging nodes. Luckily we can easily ”unhang” a hanging node by splitting the element
it lies on the edge of in two new elements, both having the hanging element in the corner, this is illustrated
in 4.6b.

E3 E4
xh xh
E1 E1 E3

E2 E2

(a) Example of a hanging node xh . (b) Unhanging of the node xh .

Figure 4.6: Example of a hanging node.

4 TWO-DIMENSIONAL WAVE EQUATION 30

4.4 Finite Element Solution

The rest of the derivation looks a lot like chapter 1 and 2. We defined the finite element solution as a linear
combination of the global basis functions like
¸
n
U px, y, tq : Ui ptqφi px, y q. (4.4.1)
i 1 
Filling this in the finite element formulation results in almost the same equations like in 1.4.6, which for
nodes xj not on the boundary are
¸
n ¸
n
: ptqT
U c2 Ui ptqSi,j  0,
i i,j

i 1
»

i 1
» (4.4.2)
with Ti,j : φi px, y qφj px, y qdA and Si,j : ∇φi px, y q  ∇φj px, y qdA.
D D

If we combine this with an equation U: ptq  0 for each node x that lies on the boundary of the domain we
j j
can construct the same ordinary differential equation as in the first two chapters, namely
: ptq
TU c2 SU ptq  0, (4.4.3)
Which we can solve in the same way as in chapter 2.

The only problem we have now is how to calculate the coefficient³ Ti,j and °n S³i,j , but luckily we can also use
here a substitution of variables. Firstly we write every integral D to i1 Ei 2 such that in each integral
over an element Ej we can write the basis functions as φi px, y q  φl,k pτj
1 px, y qq, where k is the index of the
local basis function the global basis function takes on element j, or else the basis function is zero. Then for
functions of the form f pτj
1 px, y qq, such as the integrands in Ti,j and Si,j , it holds that
» »
f pτj
1 px, y qqdA  f px, y qJdA, (4.4.4)
Ej Ec

where J is the determinant of the Jacobian matrix of τj px, y q. If the element Ej is the triangle with the
points px1 , y1 q, px2 , y2 q and px3 , y3 q at the corner, we can calculate the value J as follows

 det Bτj px, yq  Bτj px, yq  det x2
x1 x3
x1
 px2
x1 qpy3
y1 q
px3
x1 qpy2
y1 q.
Bx By y2
y1 y3
y1
J (4.4.5)

Therefore J is a constant independent on the variables of integration and can be placed outside of the integral.
³
Now the only things we need to calculate the the coefficients Ti,j and Si,j are the values Ec φl,i px, y qφl,j px, y qdA
³
and Ec ∇φl,i px, y q ∇φl,j px, y qdA (for which no actual integration is required, because the derivative are con-
stant) for i, j  1, 2, 3. These values are summarized in the tables 2 and 3.

Table 1: Values of the integrals of the local basis functions.

³ ³
Table 2: Ec
φl,i px, y qφl,j px, y qdA Table 3: Ec
∇φl,i px, y q  ∇φl,j px, y qdA
HH j HH j
H 1 2 3 H 1 2 3
i HH i HH
1 1
12
241
241 1 1
12
12
2
241 1
4
1
8 2
12 1
2 0
3
241 1
8
1
12 3
12 0 1
2

2 In cases where we can not fully tile the domain with triangles, we lose integration over the uncovered domain, but luckily

the global basis functions are defined such that they are zero on the uncovered areas.
4 TWO-DIMENSIONAL WAVE EQUATION 31

This table allows us to quickly calculate parts of the coefficients of Ti,j and Si,j , for example, if the element
Ej is the triangle with the points px1 , y1 q, px2 , y2 q and px3 , y3 q, then we get
»
φl,1 pτj
1 px, y qqφl,1 pτj
1 px, y qqdV  121 ppx2
x1 qpy3
y1 q
px3
x1 qpy2
y1 qq, (4.4.6)
Ej

and
»
∇φl,2 pτj
1 px, y qq  ∇φl,1 pτj
1 px, y qqdV  21 ppx2
x1 qpy3
y1 q
px3
x1 qpy2
y1 qq. (4.4.7)
Ej

Now that we have a way to construct matrices of the system of ordinary differential equations we need to
discretize it in time, but luckily we can use the exact same methods as we have seen in chapter 1.

4.5 Examples
We will now take a look at a few examples. First, let the constant c be 1 and the domain D be the square
r0, 2s2 . We will use the following initial values:



u0 px, y q  cos px
1q π2 cos py
1q π2 , (4.5.1)

and
u1 px, y q  0. (4.5.2)
We will split the domain regularly using nx and ny squares in the x-and y-axis respectively and divide each
square in two triangles as illustrated in figure 4.7.

0
0 2

Figure 4.7: Example discretization using nx  2 and ny  3.

In figures 4.8 we can see two examples with time step 0.01 and different amount of triangles. These figures
show for different domain discretizations the simulation at different moments in time. It is interesting to see
that after some time the shape of the wave stretches in the direction of the diagonals of the elements.
4 TWO-DIMENSIONAL WAVE EQUATION 32

(a) A 4  4 grid at time t  0. (b) A 10  10 grid at time t  0.

(c) A 4  4 grid at time t  2.5. (d) A 10  10 grid at time t  2.5.

(e) A 4  4 grid at time t  6.5. (f) A 10  10 grid at time t  6.5.

 
Figure 4.8: Simulation of the initial values u0 px, y q  cos px
1q π2 cos py
1q π2 and u1 px, y q  0 with
time step 0.01 and with pnx , ny q  p4, 4q (left) and pnx , ny q  p10, 10q (right) at different moments in time.

In an attempt to decrease the stretching we can discretize the square differently, like can be seen in figure
4.9. In this discretization we take into account the symmetric and radial properties of the initial value and
domain. The simulation behaves similarly to those in figure 4.8, but with better symmetry. There is still
some stretching to the top, which might be due to numerical errors.
4 TWO-DIMENSIONAL WAVE EQUATION 33

(a) A 4  4 grid at time t  0. (b) A 10  10 grid at time t  0.

(c) A 4  4 grid at time t  2.5. (d) A 10  10 grid at time t  2.5.

(e) A 4  4 grid at time t  6.5. (f) A 10  10 grid at time t  6.5.

 
Figure 4.9: Simulation of the initial values U0,0 px, y q  cos px
1q π2 cos py
1q π2 and U0,1 px, y q  0 on a
symmetric discretization with time step 0.01 and with pnx , ny q  p4, 4q (left) and pnx , ny q  p10, 10q (right)
at different moments in time.

We can also have more complex domains than rectangles. As the next example we will use the domain as in
figure 4.103
3 Some interesting property about this domain is explained at https://en.wikipedia.org/wiki/Hearing_the_shape_of_a_

drum.
4 TWO-DIMENSIONAL WAVE EQUATION 34

3
2

1
2

0
1 3
0 2 1 2

Figure 4.10: A more complex domain in two dimensions.

The domain has been discretized in two different ways and with initially a Gaussian peak centred around
p1.5, 0.5q. These discretizations are simulated with a time step of 0.01 and plotted at different moment in
time in figure 4.11. The wave that starts in the bottom-right corner creates smaller waves that move to the
top-left corner, just like we expect. We can also see that, just like in the previous examples, the increase of
the amount of elements slows down the propagation of the wave.
4 TWO-DIMENSIONAL WAVE EQUATION 35

(a) 28 Elements at t  0. (b) 112 Elements at t  0.

(c) 28 Elements at t  10. (d) 112 Elements at t  10.

Figure 4.11: Simulation on a non-rectangular domain with initially a Gaussian peak for two different discret-
izations with timestep 0.01.
5 CONCLUSION 36

5 Conclusion
In this thesis we have seen how to use the finite element method to numerically solve the one- and two-
dimensional wave equation where the solution is fixed at zero on the boundary of the domain and with an
initial value in time. We used the finite element method to discretize in space and a finite difference method
to discretize in time. We have seen that with this method we have a lot of freedom with both how to discretize
the domain and how to approximate the solutions.

We have also analysed the one-dimensional case both concretely and experimentally. Firstly, we have seen
that we need both the size h of the biggest element and the time steps dt to go to zero for the error to
converge to zero and that the degree of the polynomials p increases the rate of convergence of the spatial
discretizations. More specifically, that the order of convergence of the total error is O php dtq. Secondly, we
have seen that solutions of both the wave equation and the finite element solution have constant energy and
through experiments that, most of the time, the energy of the numerical solution diverges to infinity over
time, but we can reduce the initial error through refining the elements and the error over time by decreasing
the time steps. Thirdly, we have seen that there are restrictions on what kind of waves are solutions of the
wave equation and the numerical method and that for small values of h and dt the numerical restrictions
approximate those of the wave equation.

Many things we have seen could be done differently, like different types of basis functions and element shapes,
and the analysis of the one-dimensional case was also not totally rigorous everywhere as we have only really
discussed the stability of the solutions through experiments, but this is all just the tip of the iceberg with
regards to the finite element method and its variant. But, lets leave that as an adventure for the readers and
writer of this thesis.
 A IMPLEMENTATION 37

A Implementation
In this section some of the important parts of the implementions are explained.

The implementations can be found at www.github.com/jeverink/BachelorsThesis and are coded in IPy-

thon Notebooks. They can be run using either Jupyter 4 or Cocalc 5 , but Jupyter is recommended, due to
the used animation features. Each notebook is split into different cells, but it is recommended to run all cells
(Cell Ñ Run all) due to their dependencies.

A.1 One-dimensional wave equation

The notebook discussed in this section is Chapter 1.ipynb

Matrix construction
Let N and n : N 1 be the amount of elements and nodes respectively and h is the size of an element.
As we were able to calculate closed-form expressions for the coefficients in 1.4.10 for the values Ti,j and
Si,j . These expressions can be used to easily construct the matrices T and S using the following code:

Construction of T Construction of S
1 # Time coefficient matrix 1 # Space coefficient matrix
2 T = np . zeros (( n , n )); 2 S = np . zeros (( n , n ));
3 3
4 T [0 , 0] = 1; # Left boundary 4
5 T [N , N ] = 1; # Right boundary 5
6 6
7 for i in range (1 , N ): 7 for i in range (1 , N ):
8 for j in range (0 , N +1): 8 for j in range (0 , N +1):
9 if ( i == j ): 9 if ( i == j ):
10 T [i , j ] = (2.0/3.0)* h ; 10 S [i , j ] = (2.0/ h );
11 if ( abs (i - j ) == 1): 11 if ( abs (i - j ) == 1):
12 T [i , j ] = (1.0/6.0)* h ; 12 S [i , j ] = -(1.0/ h );

Both matrices start out with all zero coefficients. Firstly the first and last row are filled with the boundary
equations derived in 1.4.7 and then the rests of the coefficients are filled in as derived in 1.4.10.

A.2 Irregular elements and higher-order approximation

The notebook discussed in this section is Chapter 2.ipynb

Lagrange polynomials
As we want to be able to to approximate the solution on each element with polynomials of arbitrary degree
p we need to be able to construct the local basis functions, which are Lagrange polynomials. This can be
achieved using the following code which uses SymPy 6 for symbolic expressions:
4 www.jupyter.org/
5 www.cocalc.com/
6 www.sympy.org/
 A IMPLEMENTATION 38

Lagrange polynomials
1 x = Symbol ( ’x ’ );
2 def lagrange (p , i ):
3 width = 1.0/ p ;
4 this = i * width ;
5 num = 1; # Numerator
6 for j in range (0 , p + 1):
7 if ( j != i ):
8 num = Mul ( num , x - ( j * width ));
9 denom = 1; # Denominator
10 for j in range (0 , p + 1):
11 if ( j != i ):
12 denom = Mul ( denom , this - ( j * width ));
13 denom = Pow ( denom , -1);
14 return Mul ( num , denom );

To easier understand how this calculates the Lagrange polynomials we can rewrite its formula 2.2.4 as follow:
 
1
¹ x
xl,j ¹ ¹
φl,i pxq :   x
 xl,i
xl,j
  
¤¤
1 j p
x
1 l,i

xl,j 1
¤¤
1 j p 1
xl,j 1
¤¤
1 j p 1
. (A.2.1)

j i j i j i

This formulation is almost equivalent to the way it is calculated in the code.

Matrix construction
To construct the matrices T and S we need to calculate the coefficients Ti,j and Si,j which requires integration
over the different elements. Let h be the size of an element, phi1 and phi2 be to local basis functions and
trans and invtrans be the transformation and respectively inverse transformation of the canonical element to
the global element, then this results in the following methods:

Integration for T Integration for S

1 d e f i n t e g ( phi1 , phi2 , h ) : 1 d e f d i f f I n t e g ( phi1 , phi2 , h ,
2 p r o d u c t = Mul ( phi1 , p h i 2 ) ; 2 trans , invtrans ) :
3 r e t u r n h∗ i n t e g r a t e ( product , 3 # Transform l o c a l t o g l o b a l
4 (x , 0 , 1 ) ) ; 4 gPhi1 = p h i 1 . s u b s ( x ,
5 invtrans ) ;
6 gPhi2 = p h i 2 . s u b s ( x ,
7 invtrans ) ;
8 # Take g l o b a l d e r i v a t i v e
9 derGPhi1 = d i f f ( gPhi1 , x ) ;
10 derGPhi2 = d i f f ( gPhi2 , x ) ;
11 # Transform g l o b a l t o l o c a l
12 d e r P h i 1 = derGPhi1 . s u b s ( x ,
13 trans ) ;
14 d e r P h i 2 = derGPhi2 . s u b s ( x ,
15 trans ) ;
16 p r o d u c t = Mul ( derPhi1 , d e r P h i 2 ) ;
17 r e t u r n h∗ i n t e g r a t e ( product ,
18 (x , 0 , 1 ) ) ;

In the case of T we can directly integrate the product of the two basis functions and multiply that with the
size of the element because of 2.1.5. In the case of S we cannot directly take the derivative of the local basis
function because we require the derivative of the global basis functions and the easiest way to fix this is to
first transform the local basis functions to the global basis functions, then take their derivative with respect
 A IMPLEMENTATION 39

to x and finally transform them back to the canonical element.

Let e be a list such that for an element Ei the left external node is e[i] and the right external node is e[i+1],
then using the above defined integration methods we can now construct the matrices T and S as follows:

Construction of T Construction of S
1 # Time coefficient matrix 1 # Time coefficient matrix
2 T = np . zeros (( n , n )); 2 S = np . zeros (( nodeCount , nodeCount ));
3 3
4 for i in range (0 , N ): 4 for i in range (0 , N ):
5 h = e [ i + 1] - e [ i ]; 5 h = e [ i + 1] - e [ i ];
6 6 trans = ( h * x ) + ( e [ i ]);
7 7 invtrans = ( x - e [ i ])/ h ;
8 for j in range (0 , degree + 1): 8 for j in range (0 , degree + 1):
9 for m in range (0 , degree + 1): 9 for m in range (0 , degree + 1):
10 T [( degree * i ) + j , 10 S [( degree * i ) + j ,
11 ( degree * i ) + m ] 11 ( degree * i ) + m ]
12 += integ ( phi [ j ] , phi [ m ] , 12 += diffInteg ( phi [ j ] , phi [ m ] ,
13 h ); 13 h,
14 14 trans , invtrans );
15 15
16 for i in range (0 , n ): 16 for i in range (0 , n ):
17 T [0 , i ] = 0; 17 S [0 , i ] = 0;
18 T [n -1 , i ] = 0; 18 S [n -1 , i ] = 0;
19 T [0 , 0] = 1; 19
20 T [n -1 , n -1] = 1; 20

We start with all zero matrices and then fill it up with the coefficients Ti,j and respectively Si,j including the
first and final row. We than overwrite the the first and final row of both matrices with the respective values
for the equations of the basis functions at the boundary of the domain.

Finally we need to be able to construct the continuous numerical solution from the vector of values at the
nodes. If u is the vector at time t and pos is an x position on the domain, we can calculate the values U px, tq
using the following method:

Evaluation
1 # Evaluation
2 def ev (u , pos ):
3 for i in range (0 , N ):
4 if ( e [ i ] <= pos <= e [ i +1]):
5 h = e [ i + 1] - e [ i ];
6 localPos = ( pos - e [ i ])/ h ;
7 sum = 0;
8 for p in range (0 , degree +1):
9 sum += u [( i * degree ) + p ]* phi [ p ]. subs (x , localPos );
10 return sum ;
11 return 0;

First we search the element the x position is in and then transform it to canonical element, then by knowing
the vector and the element the position is in, we can calculate the actual value as the linear sum of the basis
functions evaluated at the local coordinate with coefficients from the vector. An alternative method would
be to simply calculate the dot product of the vector with the vector of global basis functions evaluated at
the coordinate x, but this is more expensive and a waste of calculations.
 A IMPLEMENTATION 40

A.3 Numerical analysis and experiments

The notebook discussed in this section is Chapter 3.ipynb

Energy
In chapter 3 we have seen the equation 3.2.5 which we can use to easily calculate the energy of a finite element
solution. If u and uDer are the vectors U and U9 respectively and Tc[i][j] and Sc[i][j] are the coefficients Ti,j
and Si,j respectively, then the function can be simply implemented as follows:

Energy calculation
1 # Energy function
2 def energy (u , uDer ):
3 acc = 0;
4 for i in range (0 , nodeCount ):
5 for j in range (0 , nodeCount ):
6 acc += Tc [ i ][ j ] * uDer [ i ] * uDer [ j ];
7 acc += ( c **2) * Sc [ i ][ j ] * u [ i ] * u [ j ];
8 return acc ;

The method is general enough that, if Tc and Sc are defined, the energy can be used for finite element
solutions in arbitrary dimensions.
 A IMPLEMENTATION 41

A.4 Two-dimensional wave equation

The notebook discussed in this section is Chapter 4.ipynb

Matrix construction
The biggest difference between the one-and two-dimensional case is the way the matrices T and S are
constructed. Before we construct them we first need to discuss how we store the nodes and elements. For a
node xi are x[i] and y[i] the x and y coordinates and b[i] is contains a boolean whether the node lies on the
boundary of the domain and therefore needs special treatment due to boundary conditions. Then triangles
is a list of triples containing the indices of the three corner nodes of the triangular element. Finally we let A
and Ad be the tables 2 and 3 respectively. Using these values we can implement the matrix construction as
follows:

Construction of T and S
1 T = np . zeros (( n , n ));
2 S = np . zeros (( n , n ));
3
4 for ( i1 , i2 , i3 ) in triangles :
5 inds = [ i1 , i2 , i3 ];
6 J = ( x [ i2 ] - x [ i1 ])*( y [ i3 ] - y [ i1 ]) -( x [ i3 ] - x [ i1 ])*( y [ i2 ] - y [ i1 ]);
7 for i in range (0 , 3):
8 for j in range (0 , 3):
9 T [ inds [ i ] , inds [ j ]] += J * A [ i ][ j ];
10 S [ inds [ i ] , inds [ j ]] += J * Ad [ i ][ j ];
11
12 for i in range (0 , n ):
13 if ( b [ i ]):
14 for j in range (0 , n ):
15 S [i , j ] = 0;
16 if ( i == j ):
17 T [i , j ] = 1;
18 else :
19 T [i , j ] = 0;

In the first loop, all the coefficients Ti,j and Si,j are calculated by summing the contributions of all elements.
These coefficients could then be copied to use them later for calculating the energy of a finite element solution.
The second loop is to change all rows corresponding to a node on the boundary and replace them with the
equations U9 i ptq  0.
REFERENCES I

References
[1] J. Whiteley, Finite Element Methods: A Practical Guide (Springer, 2017).
[2] S. Adjerid, Computer Methods in Applied Mechanics and Engineering 191, 4699 (2002), ISSN 0045-7825,
URL http://www.sciencedirect.com/science/article/pii/S0045782502004000.
[3] U. M. Ascher and C. Greif, A First Course in Numerical Methods (2011).

[4] M. G. Larson and F. Bengzon, The Finite Element Method: Theory, Implementation, and Applications
(Springer, 2013).
[5] L. N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations
(1996), URL http://people.maths.ox.ac.uk/trefethen/pdetext.html.

[6] W. A. Strauss, Partial Differential Equations: An Introduction (2008), 3rd ed.

[7] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods (Springer, 2008).