Introduction to spectral methods

Eric Gourgoulhon
Laboratoire de lUnivers et de ses Theories (LUTH)
CNRS / Observatoire de Paris
Meudon, France

Based on a collaboration with

Silvano Bonazzola, Philippe Grandclement, Jean-Alain Marck & Jerome Novak
Eric.Gourgoulhon@obspm.fr
http://www.luth.obspm.fr

4th EU Network Meeting, Palma de Mallorca, Sept. 2002

Plan

1. Basic principles
2. Legendre and Chebyshev expansions
3. An illustrative example
4. Spectral methods in numerical relativity

1
Basic principles

Solving a partial differential equation

Consider the PDE with boundary condition

Lu(x) = s(x),
Bu(y) = 0,

x U IRd

(1)

y U,

(2)

where L and B are linear differential operators.

Question: What is a numerical solution of (1)-(2) ?
Answer: It is a function u
which satisfies (2) and makes the residual
R := L
us
small.

What do you mean by small ?

Answer in the framework of

Method of Weighted Residuals (MWR):

Search for solutions u
in a finite-dimensional sub-space PN of some Hilbert space W
(typically a L2 space).
Expansion functions = trial functions : basis of PN : (0, . . . , N )
u
is expanded in terms of the trial functions: u
(x) =

N
X

u
n n(x)

n=0

Test functions : family of functions (0, . . . , N ) to define the smallness of the residual
R, by means of the Hilbert space scalar product:
n {0, . . . , N },

(n, R) = 0

Various numerical methods

Classification according to the trial functions n:
Finite difference: trial functions = overlapping local polynomials of low order
Finite element: trial functions = local smooth functions (polynomial of fixed degree
which are non-zero only on subdomains of U )
Spectral methods : trial functions = global smooth functions (example: Fourier series)

Various spectral methods

All spectral method: trial functions (n) = complete family (basis) of smooth global
functions
Classification according to the test functions n:
Galerkin method: test functions = trial functions: n = n and each n satisfy the
boundary condition : Bn(y) = 0
tau method: (Lanczos 1938) test functions = (most of) trial functions: n = n but
the n do not satisfy the boundary conditions; the latter are enforced by an additional
set of equations.
collocation or pseudospectral method: test functions = delta functions at special
points, called collocation points: n = (x xn).

Solving a PDE with a Galerkin method

Let us return to Equation (1).
Since n = n, the smallness condition for the residual reads, for all n {0, . . . , N },
(n, R) = 0 (n, L
u s) = 0
!

N
X

n , L
u
k k (n, s) = 0
k=0

N
X

u
k (n, Lk ) (n, s) = 0

k=0

N
X

Lnk u
k = (n, s) ,

k=0

where Lnk denotes the matrix Lnk := (n, Lk ).

Solving for the linear system (3) leads to the (N + 1) coefficients u
k of u

(3)

Solving a PDE with a tau method

Here again n = n, but the ns do not satisfy the boundary condition: Bn(y) 6= 0.
Let (gp) be an orthonormal basis of M + 1 < N + 1 functions on the boundary U and
let us expand Bn(y) upon it:
Bn(y) =

M
X

bpn gp(y)

p=0

The boundary condition (2) then becomes

Bu(y) = 0

N X
M
X

u
k bpk gp(y) = 0,

k=0 p=0

hence the M + 1 conditions:

N
X
k=0

bpk u
k = 0

0pM

Solving a PDE with a tau method (contd)

The system of linear equations for the N + 1 coefficients u
n is then taken to be the
N M first raws of the Galerkin system (3) plus the M + 1 equations above:
N
X
k=0
N
X

Lnk u
k

= (n, s)

0nN M 1

bpk u
k

= 0

0pM

k=0

The solution (
uk ) of this system gives rise to a function u
=

N
X
k=0

L
u(x) = s(x) +

M
X
p=0

p N M +p(x)

u
k k such that

10

Solving a PDE with a pseudospectral (collocation) method

This time: n(x) = (x xn), where the (xn) constitute the collocation points. The
smallness condition for the residual reads, for all n {0, . . . , N },
(n, R) = 0 ((x xn), R) = 0 R(xn) = 0 Lu(xn) = s(xn)

N
X

Lk (xn)
uk = s(xn)

(4)

k=0

The boundary condition is imposed as in the tau method. One then drops M + 1 raws
in the linear system (4) and solve the system
N
X
k=0
N
X
k=0

Lk (xn) u
k

= s(xn)

0nN M 1

bpk u
k

= 0

0pM

11

What choice for the trial functions n ?

Periodic problem : n = trigonometric polynomials (Fourier series)

Non-periodic problem : n = orthogonal polynomials

12

2
Legendre and Chebyshev expansions

13

Legendre and Chebyshev polynomials

Families of orthogonal
polynomials on [1, 1] :
Legendre polynomials:
Z 1
Pm(x)Pn(x) dx =
1

2
mn
2n + 1

Chebyshev polynomials:
Z 1
dx
Tm(x)Tn(x)
=
2
1x
1

(1 + 0n) mn
2
[from Fornberg (1998)]

3 2 1
P0(x) = 1, P1(x) = x, P2(x) = x
2
2

T0(x) = 1, T1(x) = x, T2(x) = 2x2 1

Both Legendre and Chebyshev polynomials are a subclass of Jacobi polynomials

14

Properties of Chebyshev polynomials

Definition: cos n = Tn(cos )
Recurrence relation : Tn+1(x) = 2xTn(x) Tn1(x)
Eigenfunctions of the singular Sturm-Liouville problem:

2
d p
dT
n
n
1 x2
=
Tn(x)
2
dx
dx
1x

Orthogonal family in the Hilbert space L2w [1, 1], equiped with the weight
w(x) = (1 x2)1/2 :
Z 1
(f, g) :=
f (x) g(x) w(x) dx
1

15

Polynomial interpolation of functions

Given a set of N + 1 nodes (xi)0iN in [1, 1], the Lagrangian interpolation of a
function u(x) is defined by the N -th degree polynomial:
IN u(x) =

N
X
i=0

u(xi)

N
Y
x xj
j=0
j6=i

xi xj

Cauchy theorem: there exists x0 [1, 1] such that

N
Y
1
u(x) IN u(x) =
u(N +1)(x0) (x xi)
(N + 1)!
i=0

Minimize u(x) IN u(x) independently of u minimize

N
Y
i=0

(x xi)

16

Chebyshev interpolation of functions

Note that

N
Y

(x xi) is a polynomial of degree N + 1 of the type xN +1 + aN xN +

i=0

(leading coefficient = 1).

Characterization of Chebyshev polynomials: Among all the polynomials of degree n
with leading coefficient 1, the unique polynomial which has the smallest maximum on
[1, 1] is the n-th Chebyshev polynomial divided by 2n1 : Tn(x)/2n1.
= take the nodes xi to be the N + 1 zeros of the Chebyshev polynomial TN +1(x) :
N
Y

(x xi) =

i=0

1
TN +1(x)
N
2

2i + 1

xi = cos
2(N + 1)

0iN

17

Spectral expansions : continuous (exact) coefficients

Case where the trial functions are orthogonal polynomials n in L2w [1, 1] for some
weight w(x) (e.g. Legendre (w(x) = 1) or Chebyshev (w(x) = (1 x2)1/2)
polynomials).
The spectral representation of any function u is its orthogonal projection on the space
of polynomials of degree N :
PN u(x) =

N
X

u
nn(x)

n=0

where the coefficients u

n are given by the scalar product:
1
(n, u)
u
n =
(n, n)

with (n, u) :=

The integral (5) cannot be computed exactly...

n(x) u(x) w(x) dx

1

(5)

18

Spectral expansions : discrete coefficients

The most precise way of numerically evaluating the integral (5) is given by
Gauss integration :
Z 1
N
X
f (x) w(x) dx =
wi f (xi)
1

(6)

i=0

where the xis are the N + 1 zeros of the polynomial N +1 and the coefficients wi are
Z 1
N
X
xij wj =
xi w(x) dx.
the solutions of the linear system
j=0

Formula (6) is exact for any polynomial f (x) of degree 2N + 1

Adaptation to include the boundaries of [1, 1] : x0 = 1, x1, . . . , xN 1, xN = 1
Gauss-Lobatto integration : xi = zeros of the polynomial
P = N +1 + N + N 1, with and such that P (1) = P (1) = 0. Exact for any
polynomial f (x) of degree 2N 1.

19

Spectral expansions : discrete coefficients (cont)

Define the discrete coefficients u
n to be the Gauss-Lobatto approximations of the
integrals (5) giving the u
ns :
N
X
1
wi n(xi) u(xi)
u
n :=
(n, n) i=0

(7)

The actual numerical representation of a function u is then the polynomial formed from
the discrete coefficients:
N
X
IN u(x) :=
u
nn(x) ,
n=0

instead of the orthogonal projection PN u involving the u

n .
Note: if (n) = Chebyshev polynomials, the coefficients (
un) can be computed by
means of a FFT [i.e. in N ln N operations instead of the N 2 operations of the
matrix product (7)].

20

Aliasing error
Proposition: IN u(x) is the interpolating polynomial of u through the N + 1 nodes
(xi)0iN of the Gauss-Lobatto quadrature: IN u(xi) = u(xi) 0 i N
On the contrary the orthogonal projection PN u does not necessarily pass through the
points (xi).
The difference between IN u and PN u, i.e. between the coefficients u
n and u
n, is
called the aliasing error.
It can be seen as a contamination of u
n by the high frequencies u
k with k > N , when
performing the Gauss-Lobato integration (7).

21

Illustrating the aliasing error: case of Fourier series

Alias of a sin(2x) wave by a sin(6x) wave

Alias of a sin(2x) wave by a sin(10x)

wave
[from Canuto et al. (1998)]

22

Convergence of Legendre and Chebyshev expansions

Hyp.: u sufficiently regular so that all derivatives up to some order m 1 exist.
m

Legendre:

truncation error :

interpolation error :
Chebyshev:

truncation error :

interpolation error :

C X (k)
ku kL2
kPN u ukL2 m
N
k=0
C
kPN u uk m1/2 V (u(m))
N
m
X
C
kIN u ukL2 m1/2
ku(k)kL2
N
k=0
m
C X (k)
kPN u ukL2w m
ku kL2w
N
k=0
m
C(1 + ln N ) X (k)
kPN u uk
ku k
m
N
k=0
m
C X (k)
kIN u ukL2w m
ku kL2w
N
k=0
m
X
C
kIN u uk m1/2
ku(k)kL2w
N
k=0

23

Evanescent error
From the above decay rates, we conclude that for a C function, the error in the
spectral expansion decays more rapidly than any power of 1/N . In practice, it is an
exponential decay.
Such a behavior is a key property of spectral methods and is called evanescent error.
(Remember that for a finite difference method of order k, the error decays only as
1/N k ).

24

3
An example
... at last !

25

A simple differential equation with boundary conditions

Let us consider the 1-D second-order linear (P)DE

d2u
du
x

4
+
4u
=
e
+ C,
dx2
dx

x [1, 1]

(8)

with the Dirichlet boundary conditions

u(1) = 0

and

u(1) = 0

and where C is a constant: C = 4e/(1 + e2).

The exact solution of the system (8)-(9) is
sinh 1 2x C
e +
u(x) = e
sinh 2
4
x

(9)

26

Resolution by means of a Chebyshev spectral method

Let us search for a numerical solution of (8)-(9) by means of the five first Chebyshev
polynomials: T0(x), T1(x), T2(x), T3(x) and T4(x), i.e. we adopt N = 4.
Let us first expand the source s(x) = ex + C onto the Chebyshev polynomials:
P4 s(x) =

4
X

sn Tn(x)

and

n=0

I4 s(x) =

4
X

sn Tn(x)

n=0

with
2
sn =
(1 + 0n)

dx

Tn(x)s(x)
1 x2
1

and

4
X
2
wiTn(xi)s(xi)
sn =
(1 + 0n) i=0

the xis being the 5 Gauss-Lobatto quadrature points for the

weight

1
1
w(x) = (1 x2)1/2: {xi} = { cos(i/4), 0 i 4} = 1, , 0, , 1
2
2

27

The source and its Chebyshev interpolant

2
s(x)
I4 s(x)

1.5

0.5

-0.5

-1
-1

-0.5

0.5

x
s0 = 0.03004, s1 = 1.130, s2 = 0.2715, s3 = 0.04488, s4 = 0.005474

28

Interpolation error and aliasing error

N=4 (5 Chebyshev polynomials)

0.001

P4 s(x) - s(x) (truncation error)

I4 s(x) - s(x) (interpolation error)
I4 s(x) - P4 s(x) (aliasing error)

0.0005

-0.0005

-0.001
-1

-0.5

0.5

x
sn sn = 2.0 107, 3.2 106, 4.5 105, 5.4 104, 1.0 1012

29

Matrix of the differential operator

The matrices of derivative operators with respect to the Chebyshev basis
(T0, T1, T2, T3, T4) are

d
=
dx

0
0
0
0
0

1
0
0
0
0

0
4
0
0
0

3
0
6
0
0

0
8
0
8
0

so that the matrix of the differential

is

4
0

Akl =
0
0
0

d
=

2
dx

0
0
0
0
0

0
0
0
0
0

4 0 32
0 24 0

0 0 48

0 0 0
0 0 0

d2
d
operator 2 4 + 4 Id on the r.h.s. of Eq. (8)
dx
dx

4
4
12 32
4 16 24 32

0
4
24 48

0
0
4
32
0
0
0
4

30

Resolution by means of a Galerkin method

Galerkin basis :

0(x) := T2(x) T0(x) = 2x2 2

1(x) := T3(x) T1(x) = 4x3 4x
2(x) := T4(x) T0(x) = 8x4 8x2

Each of the i satisfies the boundary conditions: i(1) = i(1) = 0. Note that the
is are not orthogonal.

1 0 1
0 1 0

0
0
Transformation matrix Chebyshev Galerkin: ki = 1

0
1
0
0
0
1
such that i(x) =

4
X

ki Tk (x).

k=0

Chebyshev coefficients and Galerkin coefficients: u(x) =

4
X
k=0

u
k Tk (x) =

2
X
i=0

i i(x)
u

The matrix ki relates the two sets of coefficients via the matrix product u
=

31

Galerkin system
For Galerkin method, the test functions are equal to the trial functions, so that the
condition of small residual writes
3
X
j = (i, s)
(i, Lu s) = 0
(i, Lj ) u
j=0

with
4 X
4
4 X
4
X
X
(i, Lj ) =
(kiTk , Llj Tl) =
kilj (Tk , LTl)
k=0 l=0

4 X
4
X
k=0 l=0

4
4 X
X
k=0 l=0

k=0 l=0

kilj (Tk ,

4
X
m=0

AmlTm) =

4 X
4
X
k=0 l=0
4

kilj

4
X

Aml(Tk , Tm)

m=0

XX

kilj (1 + 0k )Akl =
(1 + 0k )kiAkllj
2
2
k=0 l=0

32

Resolution of the Galerkin system

In the above expression appears the transpose matrix

2 0 1 0 0
h
i
Qik := t (1 + 0k )ki = 0 1 0 1 0
2 0 0 0 1
The small residual condition amounts then to solve the following linear system in

0, u
1, u
2):
u
= (u
u

QA
= Q s

4 8 8
0.331625
= 16 16
0 and Q s = 1.08544
with Q A
0
16 52
0.0655592
0 = 0.1596, u
1 = 0.09176, u
2 = 0.02949.
The solution is found to be u
:
The Chebyshev coefficients are obtained by taking the matrix product by
u
0 = 0.1891, u
1 = 0.09176, u
2 = 0.1596, u
3 = 0.09176, u
4 = 0.02949

33

Comparison with the exact solution

N=4

0.5
Exact solution
Galerkin

y = u(x)

0.4

0.3

0.2

0.1

0
-1

-0.5

0.5

sinh 1 2x
e
Exact solution: u(x) = e
e
sinh 2
1 + e2
x

34

Resolution by means of a tau method

Tau method : trial functions = test functions = Chebyshev polynomials T0, . . . , T4.
Enforce the boundary conditions by additional equations.
Since Tn(1) = (1)n and Tn(1) = 1, the boundary condition operator has the matrix

bpk =

1 1 1 1 1
1 1 1 1 1

(10)

The Tns being an orthogonal basis, the tau system is obtained by replacing the last
two rows of the matrix A by (10):

u
0
4 4
4
12 32
1
0 4 16 24 32
u

0 0
4
24 48 u
2
=
1 1
1
1
1 u
3
1 1
1
1
1
u
4

s0
s1
s2
0
0

The solution is found to be

u
0 = 0.1456, u
1 = 0.07885, u
2 = 0.1220, u
3 = 0.07885, u
4 = 0.02360.

35

Comparison with the exact solution

N=4

0.5
Exact solution
Galerkin
Tau

y = u(x)

0.4

0.3

0.2

0.1

0
-1

-0.5

0.5

sinh 1 2x
e
Exact solution: u(x) = e
e
sinh 2
1 + e2
x

36

Resolution by means of a pseudospectral method

Pseudospectral method : trial functions = Chebyshev polynomials T0, . . . , T4 and
test functions = (x xn).
The pseudospectral system is
4
X
k=0

LTk (xn) u
k = s(xn)

4 X
4
X

Alk Tl(xn) u
k = s(xn)

k=0 l=0

From a matrix point

= s, where
of view: T A u
1
1
1
1
1

1 1/ 2 0
1/ 2 1

0
1
0
1
T nl := Tl(xn) = 1

1 1/ 2
0 1/ 2 1
1
1
1
1
1

37

Pseudospectral system
To take into account the boundary conditions, replace the first row of the matrix
T A by b0k and the last row by b1k , and end up with the system

1
1
1
1
1
u
0
1
4 6.82843 15.3137 26.1421 28
u
u
4
4
0
12
12
2
4 1.17157 7.31371 2.14214
28 u
3
1
1
1
1
1
u
4

0
s(x1)
s(x2)
s(x3)
0

The solution is found to be

u
0 = 0.1875, u
1 = 0.08867, u
2 = 0.1565, u
3 = 0.08867, u
4 = 0.03104.

38

Comparison with the exact solution

N=4

0.5

y = u(x)

0.4

Exact solution
Galerkin
Tau
Pseudo-spectral

0.3

0.2

0.1

0
-1

-0.5

0.5

sinh 1 2x
e
Exact solution: u(x) = e
e
sinh 2
1 + e2
x

39

Numerical solutions with N = 6

0.5

y = u(x)

0.4

Exact solution
Galerkin
Tau
Pseudo-spectral

0.3

0.2

0.1

0
-1

-0.5

0.5

40

Exponential decay of the error with N

0
Galerkin
Tau
Pseudospectral

log10 || unum -uexact ||oo

-2
-4
-6
-8
-10
-12
-14
-16
2

10

12

14

N = number of Chebyshev polynomials -1

16

18

41

Not discussed here...

Spectral methods for 3-D problems

Time evolution
Non-linearities
Multi-domain spectral methods
Weak formulation

42

4
Spectral methods in numerical relativity

43

Spectral methods developed in Meudon

Pioneered by Silvano Bonazzola & Jean-Alain Marck (1986). Spectral methods within
spherical coordinates
1990 : 3-D wave equation
1993 : First 3-D computation of stellar collapse (Newtonian)
1994 : Accurate models of rotating stars in GR
1995 : Einstein-Maxwell solutions for magnetized stars
1996 : 3-D secular instability of rigidly rotating stars in GR

44

LORENE
Langage Objet pour la RElativite NumeriquE
A library of C++ classes devoted to multi-domain spectral methods, with adaptive
spherical coordinates.
1997 : start of Lorene
1999 : Accurate models of rapidly rotating strange quark stars
1999 : Neutron star binaries on closed circular orbits (IWM approx. to GR)
2001 : Public domain (GPL), Web page: http://www.lorene.obspm.fr
2001 : Black hole binaries on closed circular orbits (IWM approx. to GR)
2002 : 3-D wave equation with non-reflecting boundary conditions
2002 : Maclaurin-Jacobi bifurcation point in GR

45

Code for producing the figures of the above illustrative example available from Lorene
CVS server (directory Lorene/Codes/Spectral),
see http://www.lorene.obspm.fr

46

Spectral methods developed in other groups

Cornell group: Black holes

Bartnik: quasi-spherical slicing
Carsten Gundlach: apparent horizon finder
J
org Frauendiener: conformal field equations
Jena group: extremely precise models of rotating stars, cf. Marcus Ansorgs talk

47

Textbooks about spectral methods

D. Gottlieb & S.A. Orszag : Numerical analysis of spectral methods, Society for
Industrial and Applied Mathematics, Philadelphia (1977)
C. Canuto, M.Y. Hussaini, A. Quarteroni & T.A. Zang : Spectral methods in fluid
dynamics, Springer-Verlag, Berlin (1988)
B. Mercier : An introduction to the numerical analysis of spectral methods, SpringerVerlag, Berlin (1989)
C. Bernardi & Y. Maday : Approximations spectrales de problmes aux limites
elliptiques, Springer-Verlag, Paris (1992)
B. Fornberg : A practical guide to pseudospectral methods, Cambridge University
Press, Cambridge (1998)
J.P. Boyd : Chebyshev and Fourier spectral methods, 2nd edition, Dover, Mineola
(2001) [web page]