Lesson 23 Linear partial differential equations 1 We have seen - - PowerPoint PPT Presentation

Lesson 23 Linear partial differential equations

1

2

• We have seen that ODEs can be reduced to almost banded infinite-dimensional

equations using Chebyshev and ultraspherical polynomials

• These infinite-dimensional equations can be solved adaptively using Givens rota-

tions. The exact error in residual is calculated, which can be used to decide con- vergence

• This lecture we see how these ideas can be used for linear PDEs
• The resulting method works for arbitrary linear PDEs on a rectangle
• The first step is to do function approximation on a square

3

2D Function Approximation

4

• In 1D, we approximated functions using Chebyshev series:

f(x) ≈

n−1

• k=0

fkTk(x) = (T0(x), . . . , Tn−1(x))

• f0

. . . fn−1

• Functions are identified with vectors
• In 2D, we use a tensor product of Chebyshev series:

Functions are identified with matrices

5

• In 1D, we approximated functions using Chebyshev series:

f(x) ≈

n−1

• k=0

fkTk(x) = (T0(x), . . . , Tn−1(x))

• f0

. . . fn−1

• Functions are identified with vectors
• In 2D, we use a tensor product of Chebyshev series:

f(x, y) ≈

n−1

• k=0

m−1

• j=0

fkjTk(x)Tj(y) = (T0(x), . . . , Tn−1(x))

• f00

· · · f0(m−1) . . . ... . . . f(n−1)0 · · · f(n−1)(m−1)

• T0(y)

. . . Tm−1(y)

• Functions are identified with matrices

6

–1 1

• In 1D, we calculated the Chebyshev series by evaluating f(x) at the Chebyshev

points xn:

• 1.0
• 0.5

0.5 1.0

• 1.0
• 0.5

0.5 1.0

y x

• In 2D, we evaluate f(x, y) at a product of Chebyshev points xn,m := xn × x

m:

7

• In 1D, we have

Cnf(xn) ≈    ˇ f0 . . . ˇ fn1    where Cn can is viewed as an n × n matrix

• We will show that in 2D we have

8

• In 1D, we have

Cnf(xn) ≈    ˇ f0 . . . ˇ fn1    where Cn can is viewed as an n × n matrix

• We will show that in 2D we have

Cnf(xn,m)C

m ≈

   ˇ f00 · · · ˇ f0(m1) . . . ... . . . ˇ f(n1)0 · · · ˇ f(n1)(m1)   

9

• For any fixed y, we can expand

f(x, y) =

∞

• k=0

ˇ fk(y)Tk(x) where ˇ fk(y) =

• 2

π 1

−1

f(x, y)Tk(x) √ 1 − x2 x

• It follows that the smoothness of

is inherited from

• We can thus expand

10

• For any fixed y, we can expand

f(x, y) =

∞

• k=0

ˇ fk(y)Tk(x) where ˇ fk(y) =

• 2

π 1

−1

f(x, y)Tk(x) √ 1 − x2 x

• It follows that the smoothness of ˇ

fk(y) is inherited from f

• We can thus expand

ˇ fk(y) =

∞

• j=0

ˇ fkjTj(y)

11

• For any fixed y, we expand numerically

Cnf(xn, y) =:

• ˇ

f n

0 (y)

. . . ˇ f n

n1(y)

• We can thus expand
• This gives us

12

• For any fixed y, we expand numerically

Cnf(xn, y) =:

• ˇ

f n

0 (y)

. . . ˇ f n

n1(y)

• We can thus expand

Cm ˇ f n

k (xm) =:

• ˇ

f nm

k0

. . . ˇ f nm

k(m1)

• This gives us

13

• For any fixed y, we expand numerically

Cnf(xn, y) =:

• ˇ

f n

0 (y)

. . . ˇ f n

n1(y)

• We can thus expand

Cm ˇ f n

k (xm) =:

• ˇ

f nm

k0

. . . ˇ f nm

k(m1)

• This gives us

Cnf(xn, xm)C

m =

• ˇ

f n

0 (x m)

. . . ˇ f n

n1(x m)

• C

m

14

• For any fixed y, we expand numerically

Cnf(xn, y) =:

• ˇ

f n

0 (y)

. . . ˇ f n

n1(y)

• We can thus expand

Cm ˇ f n

k (xm) =:

• ˇ

f nm

k0

. . . ˇ f nm

k(m1)

• This gives us

Cnf(xn, xm)C

m =

• ˇ

f n

0 (x m)

. . . ˇ f n

n1(x m)

• C

m

= Cm ˇ f n

0 (xm), . . . , ˇ

f n

n1(xm)

• =
• ˇ

f nm

00

· · · ˇ f nm

0(m1)

. . . ... . . . ˇ f nm

(n1)0

· · · ˇ f nm

(n1)(m1)

15

Partial derivatives

16

17

18

19

20

21

22

23

• The mixed partial derivative

∂λ+µf ∂xλ∂yµ with a change of basis to C(λ) k (x) × C(µ) j

(y) series takes the form DλFD

µ

• By the same arguments, the conversion operator from

to takes the form where is the 1D conversion operator, defined in terms of that convert a series to a series

24

• The mixed partial derivative

∂λ+µf ∂xλ∂yµ with a change of basis to C(λ) k (x) × C(µ) j

(y) series takes the form DλFD

µ

• By the same arguments, the conversion operator from Tk(x)×Tj(y) to C(λ)

k (x)×

C(µ)

j

(y) takes the form S0λFS

0µ

where S0λ := Sλ1 · · · S0 is the 1D conversion operator, defined in terms of Sλ that convert a C(λ) series to a C(λ+1) series

25

26

27

0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .F. 1 …

1 4

…

• 2

3 1 6

…

• 3

8 1 8

…

1 6

• 4

15

…

1 8

• 5

24

… Å Å Å Å Å + 1

• 2

3 1 6

…

1 4

• 3

8 1 8

…

1 6

• 4

15 1 10

…

1 8

• 5

24 1 12

… Å Å Å Å Å Å Å Å Å .F. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å

=

28

Boundary conditions

29

u(−1, y) = g1(y), u(1, y) = g2(y) u(x, −1) = h1(x), u(x, 1) = h2(x) ∆u(x, y) = f(x, y) u

30

u(x, y) = (T0(x), T1(x), . . .)U

• T0(y)

T1(y) . . .

• u(−1, y)

u(1, y)

• =
• T0(−1)

T1(−1) . . . T0(1) T1(1) . . .

• U
• T0(y)

T1(y) . . .

• T (y)

31

u(x, y) = (T0(x), T1(x), . . .)U

• T0(y)

T1(y) . . .

• u(−1, y)

u(1, y)

• =
• T0(−1)

T1(−1) . . . T0(1) T1(1) . . .

• U
• T0(y)

T1(y) . . .

• =
• 1

−1 . . . 1 1 . . .

• U
• T0(y)

T1(y) . . .

• =: BU
• T0(y)

T1(y) . . .

32

u(−1, y) u(1, y)

• =

g1(y) g2(y)

• BU = G
• G = (ˇ

g1, ˇ g2) =

• ˇ

g10 ˇ g20 ˇ g11 ˇ g21 . . . . . .

33

u(x, −1) u(x, 1)

• =

h1(x) h2(x)

• UB = H
• H =

ˇ h1, ˇ h2

• =
• ˇ

h10 ˇ h20 ˇ h11 ˇ h21 . . . . . .

34

u(−1, y) = g1(y), u(1, y) = g2(y) u(x, −1) = h1(x), u(x, 1) = h2(x) ∂2u ∂x2 + ∂2u ∂y2 = f(x, y)

• BU = G,

UB = H, D2US

02 + S02UD 2

= F

35

u(−1, y) = g1(y), u(1, y) = g2(y) u(x, −1) = h1(x), u(x, 1) = h2(x) ∂2u ∂x2 + ∂2u ∂y2 + k2u = f(x, y)

• BU = G,

UB = H,

• D2 + k2S

02

• US

02 + S02UD 2

= F

Sylvester matrix equations

36

37

AX + XB = F A, B, F n × n

38

AX + XB = F A, B, F n × n AXC + DXB = F C D

• D−1AX + XBC−1 = D−1FC−1

39

AX + XB = F A, B, F n × n AXC + DXB = F C D

• D−1AX + XBC−1 = D−1FC−1

D C O(n) A B F O

• n2

40

AX + XB = F A = V ΛV −1 B = WΩW −1

−

41

AX + XB = F A = V ΛV −1 B = WΩW −1 Y = V −1XW V ΛV −1X + XWΩW −1 = V ΛY W −1 + V Y ΩW −1 = F

42

AX + XB = F A = V ΛV −1 B = WΩW −1 Y = V −1XW V ΛV −1X + XWΩW −1 = V ΛY W −1 + V Y ΩW −1 = F ΛY + Y Ω = V −1FW

43

ΛY + Y Ω =    λ1y00 · · · λ1y0(n−1) . . . ... . . . λny(n−1)0 · · · λny(n−1)(n−1)    +    ω1y00 · · · ωny0(n−1) . . . ... . . . ω1y(n−1)0 · · · ωny(n−1)(n−1)     

44

ΛY + Y Ω =    λ1y00 · · · λ1y0(n−1) . . . ... . . . λny(n−1)0 · · · λny(n−1)(n−1)    +    ω1y00 · · · ωny0(n−1) . . . ... . . . ω1y(n−1)0 · · · ωny(n−1)(n−1)    =    (λ1 + ω1)y00 · · · (λ1 + ωn)y0(n−1) . . . ... . . . (λn + ω1)y(n−1)0 · · · (λn + ωn)y(n−1)(n−1)   

45

ΛY + Y Ω =    λ1y00 · · · λ1y0(n−1) . . . ... . . . λny(n−1)0 · · · λny(n−1)(n−1)    +    ω1y00 · · · ωny0(n−1) . . . ... . . . ω1y(n−1)0 · · · ωny(n−1)(n−1)    =    (λ1 + ω1)y00 · · · (λ1 + ωn)y0(n−1) . . . ... . . . (λn + ω1)y(n−1)0 · · · (λn + ωn)y(n−1)(n−1)    ΛY + Y Ω = P yij = pij λi + ωj

46

ΛY + Y Ω =    λ1y00 · · · λ1y0(n−1) . . . ... . . . λny(n−1)0 · · · λny(n−1)(n−1)    +    ω1y00 · · · ωny0(n−1) . . . ... . . . ω1y(n−1)0 · · · ωny(n−1)(n−1)    =    (λ1 + ω1)y00 · · · (λ1 + ωn)y0(n−1) . . . ... . . . (λn + ω1)y(n−1)0 · · · (λn + ωn)y(n−1)(n−1)    ΛY + Y Ω = P yij = pij λi + ωj ΛY + Y Ω = V −1FW X = V Y W −1

47

Complexity

A = V ΛV −1 B = WΩW −1 V −1FW ΛY + Y Ω = V −1FW X = V Y W −1 O

• n3

O

• n3

O

• n3

O

• n3

O

• n3

Reducing PDEs to Sylvester equations

48

1

• 2

3 1 6

…

1 4

• 3

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= F

K 1 -1 1 -1 … 1 1 1 1 … O.U

=

U. 1 1

• 1 1

1 1

• 1 1

Å Å

= GT H

1

• 2

3 1 6

…

1 4

• 3

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= F

K 1 -1 1 -1 … 1 1 1 1 … O.U

=

U. 1 1

• 1 1

1 1

• 1 1

Å Å

= GT H Row reduction

1 2 1 1 −1 1

• 1

2 1 1 −1 1

• 1

2 1 −1 1 1

• 1

2 1 −1 1 1

1

• 2

3 1 6

…

1 4

• 3

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= F = = PT R

K 1 0 1 0 … 0 1 0 1 … O.U

U. 1 0 0 1 1 0 0 1 Å Å

1

• 2

3 1 6

…

1 4

• 3

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= F = PT

K 1 0 1 0 … 0 1 0 1 … O.U

= R

U. 1 0 0 1 1 0 0 1 Å Å

1

1 4

• 1

1 4

• D

2

D

2

1

• 2

3 1 6

…

1 4

• 3

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= F PT

K 1 0 1 0 … 0 1 0 1 … O.U

= R

U. 1 0 0 1 1 0 0 1 Å Å

1

1 4

• 1

1 4

• D

2

D

2

– –

= = R

U. 1 0 0 1 1 0 0 1 Å Å

˜ F

0 0 - 5

3

• 5

6

… 0 0

• 5

8

… 0 0

1 6

• 4

15

… 0 0

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

= = R

U. 1 0 0 1 1 0 0 1 Å Å

˜ F

0 0 - 5

3

• 5

6

… 0 0

• 5

8

… 0 0

1 6

• 4

15

… 0 0

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. 1 Å

1 4

Å

• 2

3 1 6

Å

• 3

8 1 8

Å … … … … Å

=

˜ F

0 0 - 5

3

• 5

6

… 0 0

• 5

8

… 0 0

1 6

• 4

15

… 0 0

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. Å Å

• 5

3 1 6

Å

• 5

8 1 8

Å … … … … Å

=

˜ F

0 0 - 5

3

• 5

6

… 0 0

• 5

8

… 0 0

1 6

• 4

15

… 0 0

1 8

… Å Å Å Å Å Å .U. 0 0 0 Å 0 0 0 Å 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 0 0 4 0 0 … 0 0 0 6 0 … 0 0 0 0 8 … 0 0 0 0 0 10 … Å Å Å Å Å Å Å .U. Å Å

• 5

3 1 6

Å

• 5

8 1 8

Å … … … … Å

U = U11 U12 U21 U22

• 2 × 2

2 × ∞ ∞ × 2 ∞ × ∞

=

˜ F U = U11 U12 U21 U22

• 2 × 2

2 × ∞ ∞ × 2 ∞ × ∞

• 5

3

• 5

6

…

• 5

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å .U22. 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 4 0 0 … 0 6 0 … 0 0 8 … 0 0 0 10 … Å Å Å Å Å .U22.

• 5

3 1 6

Å

• 5

8 1 8

Å … … … … Å

=

˜ F U = U11 U12 U21 U22

• 2 × 2

2 × ∞ ∞ × 2 ∞ × ∞

• 5

3

• 5

6

…

• 5

8

…

1 6

• 4

15

…

1 8

… Å Å Å Å .U22. 4 0 0 Å 0 6 0 Å 0 0 8 Å 0 0 0 10 Å … … … … Å + 4 0 0 … 0 6 0 … 0 0 8 … 0 0 0 10 … Å Å Å Å Å .U22.

• 5

3 1 6

Å

• 5

8 1 8

Å … … … … Å

Truncate to get Sylvester’s equation!

60

Recover solution

61

U22 U11, U12 U21 P =

• 1

1 · · · 1 1 · · · U11 U12 U21 U22

62

U22 U11, U12 U21 P =

• 1

1 · · · 1 1 · · · U11 U12 U21 U22

• =
• U11 +
• 1

1 · · · 1 1 · · ·

• U21

U12 +

• 1

1 · · · 1 1 · · ·

• U22

63

U22 U11, U12 U21 P =

• 1

1 · · · 1 1 · · · U11 U12 U21 U22

• =
• U11 +
• 1

1 · · · 1 1 · · ·

• U21

U12 +

• 1

1 · · · 1 1 · · ·

• U22
• P = (P1, P2) P1 2 × 2 P2 2 × ∞

U12 = P2 −

• 1

1 · · · 1 1 · · ·

• U22

64

R =

• U11

U12 U21 U22

• 1

1 1 . . . . . .

65

R =

• U11

U12 U21 U22

• 1

1 1 . . . . . .

• =
• U11 + U12
• 1

1 1 . . . . . .

• U21 + U22
• 1

1 1 . . . . . .

66

R =

• U11

U12 U21 U22

• 1

1 1 . . . . . .

• =
• U11 + U12
• 1

1 1 . . . . . .

• U21 + U22
• 1

1 1 . . . . . .

• R =
• R1

R2

• U21 = R2 − U22
• 1

1 1 . . . . . .

67

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

68

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

= R1 − U12

• 1

1 1 . . . . . .

69

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

= R1 − U12

• 1

1 1 . . . . . .

• P1 − ˜

BU21 = P1 − ˜ B

• R2 − U22 ˜

B

70

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

= R1 − U12

• 1

1 1 . . . . . .

• P1 − ˜

BU21 = P1 − ˜ B

• R2 − U22 ˜

B = R1 −

• P2 − ˜

BU22

• ˜

B = R1 − U12 ˜ B

71

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

= R1 − U12

• 1

1 1 . . . . . .

• P1 − ˜

BU21 = P1 − ˜ B

• R2 − U22 ˜

B = R1 −

• P2 − ˜

BU22

• ˜

B = R1 − U12 ˜ B ˜ BR = P ˜ B

72

U11 = P1 − 1 1 · · · 1 1 · · ·

• U21

= R1 − U12

• 1

1 1 . . . . . .

• P1 − ˜

BU21 = P1 − ˜ B

• R2 − U22 ˜

B = R1 −

• P2 − ˜

BU22

• ˜

B = R1 − U12 ˜ B ˜ BR = P ˜ B BH = GB

73

B(h1, h2) = h1(−1) h2(−1) h1(1) h2(1)

• g1

g2

• B =

g1(−1) g1(1) g2(−1) g2(1)

74

B(h1, h2) = h1(−1) h2(−1) h1(1) h2(1)

• g1

g2

• B =

g1(−1) g1(1) g2(−1) g2(1)

• h1(−1) =

x1 u(x, −1) = y1 u(−1, y) = g1(−1)

75

B(h1, h2) = h1(−1) h2(−1) h1(1) h2(1)

• g1

g2

• B =

g1(−1) g1(1) g2(−1) g2(1)

• h1(−1) =

x1 u(x, −1) = y1 u(−1, y) = g1(−1)

h2(−1) =

x1 u(x, 1) = y1 u(−1, y) = g1(1)

h1(1) =

x1 u(x, −1) = y1 u(1, y) = g2(−1)

h2(1) =

x1 u(x, 1) = y1 u(1, y) = g2(1)

76

• We warned last lecture not to use eigenvalue decomposition for non-symmetric

matrices! The algorithm for solving PDEs can be changed to use the Schur decompo- sition

• The approach only works for PDEs that are a sum of two terms

Other PDEs can be solved by forming a kronecker product matrix, with re- duces complexity to O

• n4
• Complicated geometries can be handled via domain decomposition methods (I

hope)

• Higher dimensions?