[PPT] - u + u = 0 , 2 u ( ) = u ( ) = 1 The basic idea: represent the PowerPoint Presentation

SLIDE 1

SPECTRAL METHODS

Lesson 11

SLIDE 2

2

In this lecture we consider the problem of computing solutions to linear ODEs:

a(N)(θ)u(N) + · · · + a(0)(θ)u = f(θ)

n the periodic interval

a(N)(z)u(N) + · · · + a(0)(z)u = f(z) inside the unit disk a(N)(x)u(N) + · · · + a(0)(x)u = f(x)

n the unit interval
We will impose boundary conditions so that the solution exists and is unique
For example, on the periodic interval we might solve:

u + u θ = 0, u(−π) = u(π) = 1

The basic idea: represent the infinite-dimensional differential operator with a finite dimensional

matrix

SLIDE 3

FOURIER SPECTRAL METHODS

SLIDE 4

We start with a simple equation:

u + u = f(θ), u(−π) = u(π)

Idea: represent u and f by their Fourier series:

∞

k=0

ˆ uk(ik + 1)eikθ =

∞

k=0

ˆ fkeikθ

The solution is obvious. But let’s anyways represent this in matrix form:

            ... −2i + 1 −i + 1 1 i + 1 2i + 1 ...                         . . . ˆ u−2 ˆ u−1 ˆ u0 ˆ u1 ˆ u2 . . .             =              . . . ˆ f−2 ˆ f−1 ˆ f0 ˆ f1 ˆ f2 . . .             

SLIDE 5

Numerically, we truncate the operator and use the DFT to compute the coefficients:
α + 1

... β + 1

um = Ff m

where f m = f(θm) =

f(θ1)

. . . f(θm)

, and
um ≈
ˆ

uα . . . ˆ uβ

In this simple case, convergence follows from convergence of fm

SLIDE 6

An alternative is to represent the operator as acting on function values:
Which method is better? Let’s compare

F−1    α + 1 ... β + 1    Fum = f m where um ≈ u(θm)

SLIDE 7

coefficient space value space

Matrix density

O(m) operations O(m3) operations

SLIDE 8

We introduce some notation
The differentiation operator is the bi-infinite operator

D =

...

−

...
The truncation operator maps bi-infinite vectors to finite vectors:

Pm u =

uα

. . . uβ

The truncated differentiation matrix is:

Dm = PmDP

m =

α

... β

.

SLIDE 9

We want to now incorporate variable coefficients, i.e., equations of the form

u + a(θ)u = f(θ)

In value space this is:
F1DmF +
a(θ1)

... a(θm)

um = f m

where Dm is the discretization of the differentiation operator: Dm =

α

... β

.
In coefficient space we have:
Dm + F
a(θ1)

... a(θm)

F1
um = Ff m

SLIDE 10

50 100 150 200 250 300 350 10-15 10-12 10-9 10-6 0.001 1

Error Number of unknowns coefficient space value space

u(θ) + u(θ) θ = 1 u(−π) = u(π)

SLIDE 11

coefficient space value space

Matrix density

O(C m) operations O(m3) operations

SLIDE 12

MULTIPLICATION OF FOURIER SERIES

SLIDE 13

We want to simplify

F

a(θ1)

... a(θm)

F−1,

i.e., we want to investigate multiplication in coefficient space

Let's consider a(θ) = −θ: we have

−θ

∞

k=−∞

ˆ fkkθ =

∞

k=−∞

ˆ fk+1kθ

Thus in operator space we have the shift operator:

... ... so that . . . . . . . . . . . .

SLIDE 14

We want to simplify

F

a(θ1)

... a(θm)

F−1,

i.e., we want to investigate multiplication in coefficient space

Let's consider a(θ) = −θ: we have

−θ

∞

k=−∞

ˆ fkkθ =

∞

k=−∞

ˆ fk+1kθ

Thus in operator space we have the shift operator:

S =

...

1 1 ...

, so that S
.

. . ˆ f−1 ˆ f0 ˆ f1 . . .

=
.

. . ˆ f0 ˆ f1 ˆ f2 . . .

SLIDE 15

More generally, if a(θ) = jθ: we have

−jθ

∞

k=−∞

ˆ fkkθ =

∞

k=−∞

ˆ fk+jkθ

This is equivalent to

times

, i.e.,

... ...

SLIDE 16

j times

More generally, if a(θ) = jθ: we have

−jθ

∞

k=−∞

ˆ fkkθ =

∞

k=−∞

ˆ fk+jkθ

This is equivalent to

j times

−θ · · · −θ, i.e.,

Sj =

...

1 · · · 1 ...

SLIDE 17

For general, smooth and periodic a we can expand in Fourier series:

a(θ) =

∞

k=−∞

ˆ akkθ

Each term is multiplication by jθ
Thus, in operator space we have the multiplication operator

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

A bi-infinite matrix with constant diagonals is called a Laurent operator

SLIDE 18

For general, smooth and periodic a we can expand in Fourier series:

a(θ) =

∞

k=−∞

ˆ akkθ

Each term is multiplication by jθ
Thus, in operator space we have the multiplication operator

L[a] =

∞

k=−∞

ˆ akS−k =

...

... ... ... ... ... ˆ a0 ˆ a−1 ˆ a−2 ... ... ˆ a1 ˆ a0 ˆ a−1 ... ... ˆ a2 ˆ a1 ˆ a0 ... ... ... ... ... ...

A bi-infinite matrix with constant diagonals is called a Laurent operator

SLIDE 19

In practice, we can approximate a:

a(θ) ≈

β

k=α

ˆ am

k kθ

Thus, in operator space the multiplication operator is banded
For example, if

we have

SLIDE 20

In practice, we can approximate a:

a(θ) ≈

β

k=α

ˆ am

k kθ

Thus, in operator space the multiplication operator is banded
For example, if α = β = 1 we have

L[a] =

...

... ... ... ... ... ˆ a0 ˆ a−1 ˆ a−2 ... ... ˆ a1 ˆ a0 ˆ a−1 ... ... ˆ a2 ˆ a1 ˆ a0 ... ... ... ... ... ...

≈
...

... ... ˆ am ˆ am

−1

ˆ am

1

ˆ am ˆ am

−1

ˆ am

1

ˆ am ... ... ...

SLIDE 21

Finally, we truncate the multiplication operator. Thus the operator

u + a(θ)u = f becomes (Dm + Lm[a]) um = Ff m where Lm[a] = PmL[a]P

m =

ˆ

a0 ˆ a1 ˆ a1 ˆ a0 ˆ a1 ... ... ... ˆ a1 ˆ a0 ˆ a1 ˆ a1 ˆ a0

SLIDE 22

HIGHER ORDER EQUATIONS

SLIDE 23

The operator

a(N)(θ)u(N) + · · · + a(0)u in coefficient space is L[a(N)]DN + · · · L[a(0)]

We can truncate it by

Lm[a(N)]DN

m + · · · Lm[a(0)]

However, we need to impose boundary conditions so that the operator is uniquely

invertible

In finite dimensions, non-unique inverse causes bad conditioning

SLIDE 24

We denote the boundary conditions as an d × ∞ operator:

B u

An example: suppose we want to impose Neumann conditions at the point χ:

u(χ) = u(χ) = 0 Then we have the boundary condition operator B =

· · ·

2χ χ 1 χ 2χ · · · · · · −22χ −χ χ 22χ · · ·

SLIDE 25

The differential equation

Bu = c and a(N)(θ)u(N) + · · · + a(0)(θ)u = f(θ) in coefficient space is a pair of equations: B u = c and

L[a(N)]DN + · · · + L[a(0)]
u =

f

We need a square matrix, so we now truncate it as
BP

m

Pmd

L[a(N)]DN + · · · + L[a(0)]
P

m

um =
c

Ff md

SLIDE 26

TAYLOR SERIES SPECTRAL METHODS

SLIDE 27

For Taylor series, the vector of coefficients becomes only infinite in one direction
u =
ˆ

u0 ˆ u1 . . .

Differentiation becomes

D =

1

2 3 ...

Truncation is now defined by

Pm u =

ˆ

u0 . . . ˆ um−1

SLIDE 28

j times

Multiplication by zj is again a shift operator:

T[zj] =           . . . 1 1 ...          

Thus multiplication by a(z) becomes a Toeplitz operator

T[a] =      ˆ a0 ˆ a1 ˆ a0 ˆ a2 ˆ a1 ˆ a0 . . . ... ... ...     

SLIDE 29

The differential equation

a(N)(z)u(N) + · · · + a(0)(z)u = f(z) and Bu = c becomes

BP

m

Pmd

T[a(N)]DN + · · · + T[a(0)]
P

m

um =
c

T f md

SLIDE 30

Example: the Airy function (z) satisfies

(z) = z(z) and (0) = 1 32/3Γ 2

3

, (0) = 1 31/3Γ 1

3

The operator becomes

... ...

We truncate the operator and solve the associated equation:

... ...

.

. .

SLIDE 31

Example: the Airy function (z) satisfies

(z) = z(z) and (0) = 1 32/3Γ 2

3

, (0) = 1 31/3Γ 1

3

The operator becomes
B

D2 + T[z]

=
1

1 2 1 6 1 12 ... ...

We truncate the operator and solve the associated equation:

... ...

.

. .

SLIDE 32

Example: the Airy function (z) satisfies

(z) = z(z) and (0) = 1 32/3Γ 2

3

, (0) = 1 31/3Γ 1

3

The operator becomes
B

D2 + T[z]

=
1

1 2 1 6 1 12 ... ...

We truncate the operator and solve the associated equation:
1

1 2 1 6 ... ... 1 (m − 1)!

um =
(0)

(0) . . .

SLIDE 33

Convergence at 1

10 20 30 40 50 60 10-14 10-11 10-8 10-5 0.01 10 20 30 40 50 60 10-14 10-11 10-8 10-5 0.01

Convergence at 2

10 20 30 40 50 60 10-44 10-35 10-26 10-17 10-8

Computed Taylor coefficients

SLIDE 34

CHEBYSHEV SPECTRAL METHODS

SLIDE 35

Recall from last lecture we found a fast way of computing a derivative of function
In other words, we can apply the differentiation operator Dˇ

u efficiently

On the other hand, the differentiation operator itself is dense:

D =              1 3 5 · · · 4 8 12 6 10 ... 8 12 10 ... 12 ...             

Thus imposing the equation in coefficient space loses efficiency, and hence we will

focus on constructing the equation in value space Next lecture we will adapt the approach to regain banded matrices

SLIDE 36

Recall the discrete cosine transform C which maps function values at Chebyshev

points f n = f(xn) to (approximate) Chebyshev coefficients Cf n ≈    ˇ f0 . . . ˇ fn−1   

Using this, we can represent differentiation in value space by

Dn = C−1DnC

While this can be slow to generate for each n, the computational can be reused

for multiple ODEs

SLIDE 37

The second-order differential operator

a(2)(x)u + a(1)(x)u + a(0)(x)u becomes

(a(2)(xn))D2

n + (a(1)(xn))Dn + (a(0)(xn))

un
There's a problem: which rows do we drop to incorporate the boundary condi-

tions?

The typical boundary conditions are at the endpoints ±1

For example, Dirichlet conditions are Bun =

e

1

e

1

un =
1

· · · · · · 1

un =
u(x1)

u(xn)

Thus the standard approach is to drop the boundary rows

SLIDE 38

The second-order differential equation

a(2)(x)u + a(1)(x)u + a(0)(x)u = f(x) and Bu = c becomes

BP

n

I2:n1

(a(2)(xn))D2

n + (a(1)(xn))Dn + (a(0)(xn))

un =
c

f n

for the n − 2 × n boundary row dropping matrix

I2:n1 =

1

... 1