In other words, we can apply the differentiation operator D u - - PowerPoint PPT Presentation

SLIDE 1

ULTRASPHERICAL SPECTRAL METHODS

Lesson 12

SLIDE 2

• 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 We will adapt the approach to regain banded matrices

SLIDE 3

VALUE-SPACE CHEBYSHEV SPECTRAL METHODS

SLIDE 4

• 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 5

• 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 6

• 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

SLIDE 7

Value space Chebyshev spectral method

50 100 150 200 10-12 10-10 10-8 10-6 10-4 0.01

Error

u00(x) + exu(x) = 0, u(−1) = 0, u(1) = 1

Matrix sparsity

SLIDE 8

           

• −2

−

• 2
• 

          

         1 2 3 4 5

• 

       

             1 3 5 · · · 4 8 12 6 10

• 8

12 10

• 12
• 

           

L[a] =              

• ˆ

a0 ˆ a−1 ˆ a−2 ˆ a−3

• ˆ

a1 ˆ a0 ˆ a−1 ˆ a−2

• ˆ

a2 ˆ a1 ˆ a0 ˆ a−1

• ˆ

a3 ˆ a2 ˆ a1 ˆ a0

• 

            

T[a] =        ˆ a0 ˆ a1 ˆ a0 ˆ a2 ˆ a1 ˆ a0 ˆ a3 ˆ a2 ˆ a1 ˆ a0

• 

     

Fourier Taylor Chebyshev

Differentiation Multiplication

?

SLIDE 9

Fourier Taylor Chebyshev

Differentiation Multiplication

SLIDE 10

FIRST ORDER DIFFERENTIATION

SLIDE 11

• We want to modify our method for Chebyshev polynomials so that differentiation is

banded

• The idea: represent the derivative in a different basis
• We use the Chebyshev U polynomials

– We define these by differentiating the Chebyshev T polynomials:

Uk(x) = T

k+1(x)

k + 1 , for k = 0, 1, 2, . . .

SLIDE 12

• Thus if

a(x) =

• k=0

ˇ akTk(x) we represent the derivative as a(x) =

• k=0

ˇ ak+1T

k+1(x) =

• k=0

ˇ ak+1(k + 1)Uk(x)

• In matrix form we get the differentiation operator from T to U

D =

• 1

2 3 4 ...

SLIDE 13

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) =

• x
• Tk(x) x

SLIDE 14

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) =

• x
• Tk(x) x

=

• x

Tk+1(x) 2(k + 1) − Tk1(x) 2(k − 1) + C

• = T

k+1(x)

2(k + 1) − T

k1(x)

2(k − 1)

SLIDE 15

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) =

• x
• Tk(x) x

=

• x

Tk+1(x) 2(k + 1) − Tk1(x) 2(k − 1) + C

• = T

k+1(x)

2(k + 1) − T

k1(x)

2(k − 1) = Uk(x) 2 − Uk2(x) 2

SLIDE 16

• We have the conversion relationships

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) = Uk(x) 2 − Uk−2(x) 2

• We use these to construct the T to U conversion operator:
• In other words, if

and then

SLIDE 17

• We have the conversion relationships

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) = Uk(x) 2 − Uk−2(x) 2

• We use these to construct the T to U conversion operator:

S0 =

• 1

− 1

2 1 2

− 1

2 1 2

− 1

2

... ...

• In other words, if

and then

SLIDE 18

• We have the conversion relationships

T0(x) = U0(x), T1(x) = U1(x) 2 , Tk(x) = Uk(x) 2 − Uk−2(x) 2

• We use these to construct the T to U conversion operator:

S0 =

• 1

− 1

2 1 2

− 1

2 1 2

− 1

2

... ...

• In other words, if

f(x) =

∞

• k=0

ˇ fkTk(x) and

• u0

u1 . . .

• = U
• ˇ

f0 ˇ f1 . . .

• then

f(x) =

∞

• k=0

ukUk(x)

SLIDE 19

• Example: u + u = f and u(−1) = 1
• We represent u in Chebyshev T series and f in Chebyshev U series
• Thus the ODE becomes

(D + S0)ˇ u = S0ˇ f

• The boundary conditions are
• Thus we define the functional

for

SLIDE 20

• Example: u + u = f and u(−1) = 1
• We represent u in Chebyshev T series and f in Chebyshev U series
• Thus the ODE becomes

(D + S0)ˇ u = S0ˇ f

• The boundary conditions are

u(−1) =

• k=0

ˇ ukTk(−1) =

• k=0

ˇ uk k (−1)

• Thus we define the functional

for

SLIDE 21

• Example: u + u = f and u(−1) = 1
• We represent u in Chebyshev T series and f in Chebyshev U series
• Thus the ODE becomes

(D + S0)ˇ u = S0ˇ f

• The boundary conditions are

u(−1) =

• k=0

ˇ ukTk(−1) =

• k=0

ˇ uk k (−1) =

• k=0

ˇ uk kπ =

• k=0

ˇ uk(−1)k

• Thus we define the functional

for

SLIDE 22

• Example: u + u = f and u(−1) = 1
• We represent u in Chebyshev T series and f in Chebyshev U series
• Thus the ODE becomes

(D + S0)ˇ u = S0ˇ f

• The boundary conditions are

u(−1) =

• k=0

ˇ ukTk(−1) =

• k=0

ˇ uk k (−1) =

• k=0

ˇ uk kπ =

• k=0

ˇ uk(−1)k

• Thus we define the functional

Bˇ u = 1 for B = (1, −1, 1, −1, . . .)

SLIDE 23

• Example: u + u = f and u(−1) = 1
• This becomes
• B

D + S0

• ˇ

u =

• 1

−1 1 −1 1 · · · 1 1 − 1

2 1 2

2 − 1

2 1 2

3 − 1

2

... ... ...

• ˇ

u = 1 S0ˇ f

SLIDE 24

MULTIPLICATION OF CHEBYSHEV SERIES

SLIDE 25

L[a] =              

• ˆ

a0 ˆ a−1 ˆ a−2 ˆ a−3

• ˆ

a1 ˆ a0 ˆ a−1 ˆ a−2

• ˆ

a2 ˆ a1 ˆ a0 ˆ a−1

• ˆ

a3 ˆ a2 ˆ a1 ˆ a0

• 

            

T[a] =        ˆ a0 ˆ a1 ˆ a0 ˆ a2 ˆ a1 ˆ a0 ˆ a3 ˆ a2 ˆ a1 ˆ a0

• 

     

Fourier Taylor Chebyshev

Relationship Operator

?

zjzk = zj+k

jθkθ = (j+k)θ

Tj(x)Tk(x) = ?

SLIDE 26

• We go to the circle:

Tj(J(z))Tk(J(z)) = zj + z−j 2 zk + z−k 2

• In other words,

SLIDE 27

• We go to the circle:

Tj(J(z))Tk(J(z)) = zj + z−j 2 zk + z−k 2 = zj+k + zj−k + zk−j + z−j−k 4

• In other words,

SLIDE 28

• We go to the circle:

Tj(J(z))Tk(J(z)) = zj + z−j 2 zk + z−k 2 = zj+k + zj−k + zk−j + z−j−k 4 = zj+k + z−j−k 4 + zj−k + zk−j 4

• In other words,

SLIDE 29

• We go to the circle:

Tj(J(z))Tk(J(z)) = zj + z−j 2 zk + z−k 2 = zj+k + zj−k + zk−j + z−j−k 4 = zj+k + z−j−k 4 + zj−k + zk−j 4 = 1 2 zj+k + z−j−k 2 + 1 2 z|j−k| + z−|j−k| 2

• In other words,

SLIDE 30

• We go to the circle:

Tj(J(z))Tk(J(z)) = zj + z−j 2 zk + z−k 2 = zj+k + zj−k + zk−j + z−j−k 4 = zj+k + z−j−k 4 + zj−k + zk−j 4 = 1 2 zj+k + z−j−k 2 + 1 2 z|j−k| + z−|j−k| 2 = Tj+k(J(z)) 2 + T|j−k|(J(z)) 2

• In other words,

Tj(x)Tk(x) = Tj+k(x) 2 + T|j−k|(x) 2

SLIDE 31

• We now construct the multiplication operators in coefficient space term by term
• Since T0(x) = 1, the first case is trivial:

M[T0] =    1 1 ...   

• Now for

we have ... ... . . . . . . ...

SLIDE 32

• We now construct the multiplication operators in coefficient space term by term
• Since T0(x) = 1, the first case is trivial:

M[T0] =    1 1 ...   

• Now for T1(x) we have

M[T1] = 1 2      1 1 1 1 1 ... ...      + 1 2      · · · 1 · · · · · · . . . . . . ...     

SLIDE 33

• M[T2] = 1

2        1 1 1 1 1 1 ... ...        + 1 2       1 1      

• ...

...

• We now see a pattern

SLIDE 34

• M[T2] = 1

2        1 1 1 1 1 1 ... ...        + 1 2       1 1      

• M[T3] = 1

2          1 1 1 1 1 1 1 ... ...          + 1 2       1 1 1      

• We now see a pattern

SLIDE 35

M[a] = 1 2                     2a0 a1 a2 a3 . . . a1 2a0 a1 a2 ... a2 a1 2a0 a1 ... a3 a2 a1 2a0 ... . . . ... ... ... ...           +          . . . a1 a2 a3 a4 . . . a2 a3 a4 a5 ... a3 a4 a5 a6 ... . . . ... ... ... ...                    .

Multiplication operator (Toeplitz + Hankel)

a(x) =

∞

• k=0

akTk(x)

SLIDE 36

Multiplication operator (Toeplitz + Hankel) Because approximation of a by its Chebyshev series converges spectrally fast, this matrix is in practice banded

M[a] = 1 2                   2a0 a1 a2 a1 2a0 a1 a2 a2 a1 2a0 a1 ... a2 a1 2a0 ... ... ... ...          +         . . . a1 a2 a2                  . a(x) =

m

• k=0

akTk(x)

SLIDE 37

• Our operator equation representing

u + a(x)u = f and u(−1) = 1 is now

• B

D + S0M[a]

• ˇ

u = 1 S0ˇ f

• This is banded

SLIDE 38

Computed Solution

u + x3u = 100 sin(20,000x2), u(−1) = 0

Error at one

5000 10000 15000 20000 25000 10-10 10-8 10-6 10-4 0.01 1

SLIDE 39

u + sin(2,000x2)u = 0, u(−1) = 1

−1 −0.5 0.5 1 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005 1.01

x u(x)

Solution

2000 4000 6000 8000 10000 12000 10

−30

10

−25

10

−20

10

−15

10

−10

10

−5

10 ||un−un−1||2 n

Cauchy Error

SLIDE 40

SECOND ORDER ODES

SLIDE 41

• We replicate the approach of first order equations:

– We construct a new basis so that second order differentiation is banded – We find a conversion relationship between the Chebyshev T basis and the new basis – We use the bandedness of multiplication in coefficient space

SLIDE 42

• Step 1: define a new basis, the ultraspherical polynomials

C(2)

k (x) = U k+1(x)

2

• Then second order differentiation is

D2 = 2      1 1 1 ...      D = 2      1 1 1 ...             1 2 3 4 ...        = 2      2 3 4 ...     

SLIDE 43

• Step 2: We need to convert Chebyshev U polynomials to ultraspherical C poly-

nomials

• The first two can be found directly:

U0(x) = C(2)

0 (x), U1(x) = C(2) 1 (x)

2

• We now use the same trick of integrating then differentiating:

Uk(x) =

• x
• Uk(x) x =
• x

Tk+1(x) k + 1

SLIDE 44

• Step 2: We need to convert Chebyshev U polynomials to ultraspherical C poly-

nomials

• The first two can be found directly:

U0(x) = C(2)

0 (x), U1(x) = C(2) 1 (x)

2

• We now use the same trick of integrating then differentiating:

Uk(x) =

• x
• Uk(x) x =
• x

Tk+1(x) k + 1 = 1 k + 1

• x

Uk+1(x) 2 − Uk−1(x) 2

SLIDE 45

• Step 2: We need to convert Chebyshev U polynomials to ultraspherical C poly-

nomials

• The first two can be found directly:

U0(x) = C(2)

0 (x), U1(x) = C(2) 1 (x)

2

• We now use the same trick of integrating then differentiating:

Uk(x) =

• x
• Uk(x) x =
• x

Tk+1(x) k + 1 = 1 k + 1

• x

Uk+1(x) 2 − Uk−1(x) 2

• = C(2)

k (x) − C(2) k−2(x)

k + 1

SLIDE 46

• We have the conversion relationships:

U0(x) = C(2)

0 (x), U1(x) = C(2) 1 (x)

2 , Uk(x) = C(2)

k (x) − C(2) k−2(x)

k + 1

• Thus we can construct the banded operator that converts Chebyshev U to Ultra-

spherical C:

SLIDE 47

• We have the conversion relationships:

U0(x) = C(2)

0 (x), U1(x) = C(2) 1 (x)

2 , Uk(x) = C(2)

k (x) − C(2) k−2(x)

k + 1

• Thus we can construct the banded operator that converts Chebyshev U to Ultra-

spherical C: S1 =      1 − 1

3 1 2

− 1

4 1 3

− 1

5

... ...     

SLIDE 48

• We want to solve the ODE

u + a(x)u = f(x), u(−1) = c1 and u(1) = c1

• We again represent the unknown u in Chebyshev T series but now the right-hand

side in ultraspherical C series

• For the multiplication operator, we multiply in T series, then convert to U series,

then to C series

• Thus we get
• B

D2 + S1S0M[a]

• ˇ

u =

• c1

c1 S1S0ˇ f

SLIDE 49

D2

+

M[−x] u = xu, u(−1) = Ai (−1), u(1) = Ai (1)

Example:

S1 S0

SLIDE 50

1 0 - 1

3 1 2

• 1

4

0 0

1 3

• 1

5

0 0

1 4

• 1

6

0 0

1 5

• 1

7

0 0

1 6

• 1

8

0 0

1 7

• 1

9

0 0

1 8

• 1

10

0 0

1 9

0 0

1 10

1 0 - 1

2 1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

• 1

2

0 0

1 2

0 0

1 2

• 1

2

• 1
• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

• 1

2

0 0 4 0 0 0 0 0 6 0 0 0 0 0 8 0 0 0 0 0 10 0 0 0 0 0 12 0 0 0 0 0 14 0 0 0 0 0 16 0 0 0 0 0 18 0 0 0 0 0 0 0 0 0 0

D2

+

M[−x] u = xu, u(−1) = Ai (−1), u(1) = Ai (1)

Example:

S1 S0

SLIDE 51

1

• 1

1

• 1

1

• 1

1

• 1

1

• 1

1 1 1 1 1 1 1 1 1 1

• 1

6

4

1 4

• 1

12

• 1

4 1 16

6

1 8

• 1

16

• 1

12 1 20

8

1 12

• 1

20

• 1

16 1 24

10

1 16

• 1

24

• 1

20 1 28

12

1 20

• 1

28

• 1

24 1 32

14

1 24

• 1

28 1 36

16

4 63

• 1

32 1 40

18

• B

D2 + S1S0M[−x]

• =

SLIDE 52

20 40 60 80 100 10-112 10-89 10-66 10-43 10-20

u = xu, u(−1) = Ai (−1), u(1) = Ai (1) Computed ˇ u100

Computed coefficients

SLIDE 53

Convergence of derivatives

20 40 60 80 100 120 140 10-13 10-10 10-7 10-4 0.1

u = xu, u(−1) = Ai (−1), u(1) = Ai (1)

n d = 0, 9, 20, 99

• u(d)

n (1/2) − Ai(d)(1/2)

SLIDE 54

• 1.0
• 0.5

0.5 1.0

• 8
• 6
• 4
• 2

2 4

u + 5002 ex3 − 1

• u = sin(500x2), u(−1) = 1, u(1) = −1

Solution Cauchy error

100 200 300 400 500 600 700 10-11 10-8 10-5 0.01 10