Laurette TUCKERMAN laurette@pmmh.espci.fr Numerical Methods for - - PowerPoint PPT Presentation

Laurette TUCKERMAN laurette@pmmh.espci.fr

Numerical Methods for Differential Equations in Physics

Fast Fourier Transform

ˆ uk+N =

N−1

• ℓ=0

u(xℓ)e−i(k+N)2πℓ/N =

• u(xℓ)e−ik2πℓ/Ne−iN2πℓ/N

=

• u(xℓ)e−ik2πℓ/N = ˆ

uk ˆ uk is N-periodic ˆ uk =

N/2−1

• ℓ=0

u(x2ℓ)e−i2π(2ℓ)k/N +

N/2−1

• ℓ=0

u(x2ℓ+1)e−i2π(2ℓ+1)k/N =

N/2−1

• ℓ=0

u(x2ℓ)e−i2πℓk/(N/2) + e−i2πk/N

N/2−1

• ℓ=0

u(x2ℓ+1)e−i2πℓk/(N/2) =

N/2−1

• ℓ=0

v(xℓ)e−i2πℓk/(N/2) + e−i2πk/N

N/2−1

• ℓ=0

w(xℓ)e−i2πℓk/(N/2) = ˆ vk + e−i2πk/N ˆ wk ˆ uk+N/2 = ˆ vk − e−i2πk/N ˆ wk v, w are of length N/2

Two transforms each of size N/2, and N additions → 2(N/2)2 + N Each trans of size N/2 becomes two trans of size N/4, and N/2 adds. Each trans of size N/4 becomes two trans of size N/8, and N/4 adds. N 2 − → 2 N 2 2 + N = N 2 2 + N N 4 − → 2

• 2

N 4 2 + N 2

• + N = N 2

4 + 2N N 8 − → 2

• 2
• 2

N 8 2 + N 4

• + N

2

• + N = N 2

8 + 3N N = 2p: N 2p − → N 2 2p + pN = N + (log2 N)N = O(N log2 N)

Fourier transform in one dimension:

• uk =

N−1

• ℓ=0

u(xℓ)e−i2πkℓ/N One multiplication for each k, ℓ, so O(N 2

x) operations.

Fourier transform in two dimensions:

• uk,m =

Ny−1

• n=0

Nx−1

• ℓ=0

u(xℓ, yn)e−i2πkℓ/Nx

• uk(yn)

e−i2πmn/Ny {u(xℓ, yn)}

N 2

xNy

− − − → { uk(yn)}

NxN 2

y

− − − → {

• uk,m}

Total: NxNy(Nx + Ny) ≪ N 2

xN 2 y

With FFT: NxNy(log Nx + log Ny) = NxNy log(NxNy)

Even without FFT, multidimensional Fourier transform would be fast be- cause the different dimensions are decoupled: NxNyNz(Nx + Ny + Nz) Fourier transform in x is action with matrix F x

k,ℓδm,nδi,j

Fourier transform in y is action with matrix F y

m,nδk,ℓδi,j

Fourier transform in z is action with matrix F z

i,jδm,nδk,ℓ

Contrast with convolution: not separable Convolution in one dimension

• fg(k) =
• f ∗

g

• (k) =
• ℓ
• fℓ

gk−ℓ One multiplication for each k, ℓ, so O(N 2

x) operations.

Convolution in two dimensions:

• fg(k, m) =
• f ∗

g

• (k, m) =
• ℓ
• n
• fℓ,n

gk−ℓ,m−n One multiplication for each k, ℓ, m, n, so O(N 2

xN 2 y) operations.

Multidimensional Fourier transform of real data

For u(x) real, ˆ uk is conjugate symmetric: ˆ u−k = ˆ u∗

k so that

u(x) ∼ ˆ ukeikx + ˆ u−ke−ikx = ˆ ukeikx +

• ˆ

ukeikx∗ = 2Re

• ˆ

ukeikx For u(x, y) real, ˆ u−k,−m = ˆ u∗

k,m so need half of (k, m) plane.

• r

For u(x, y, z) real, ˆ u−k,−m,−n = ˆ u∗

k,m,n so need half of (k, m, n) plane.

Discretization in Two or Three Dimensions

Periodic BCs in x and y: Fourier-Fourier

u(x, y) =

• k,m

ˆ uk,meikxeimy g(x, y) =

• k,m

ˆ gk,meikxeimy ∆u =

• k,m

(−k2 − m2)ˆ uk,meikxeimy =

• k,m

ˆ gk,meikxeimy ˆ uk,m = −ˆ gk,m k2 + m2 Have used interval [0, 2π) for simplicity. More generally, domain is [0, Lx) × [0, Ly) and basis functions are eikx2π/Lxeimy2π/Ly

Periodic BCs in x, Dirichlet BCs in y: Fourier-Finite Differences

u(x, y) =

• k

ˆ uk(y)eikx g(x, y) =

• k

ˆ gk(y)eikx ∆u =

• k
• −k2ˆ

uk(y) + ˆ uk(y + ∆y) − 2ˆ uk(y) + ˆ uk(y − ∆y) (∆y)2

• eikx

=

• k

ˆ uk(y + ∆y) +

• −2 − (k∆y)2

ˆ uk(y) + ˆ uk(y − ∆y) (∆y)2 eikx For each Fourier component k

1 (∆y)2        −2 − (k∆y)2 1 1 −2 − (k∆y)2 1 1 −2 − (k∆y)2 1 ... 1 −2 − (k∆y)2)               ˆ uk(y1) ˆ uk(y2) ˆ uk(y3) . . . ˆ uk(yN)       

Boundary conditions needed for each Fourier component k, e.g.

1 (∆y)2        1 . . . 1 −2 − (k∆y)2 1 1 −2 − (k∆y)2 1 ... . . . 1               ˆ uk(y1) ˆ uk(y2) ˆ uk(y3) . . . ˆ uk(yN)        =        γ− ˆ gk(y2) ˆ gk(y3) . . . γ+       

Elliptic equations need boundary conditions (Dirichlet, Neumann or peri-

• dic) along all boundaries of domain

Periodic BCs in x, Dirichlet BCs in y: Fourier-Chebyshev

u(x, y) =

• k,n

uk,neikxTn(y) g(x, y) =

• k,n

gk,neikxTn(y) ∆u =

• k,n

−k2uk,neikxTn(y) +

• k,n

uk,neikxT ′′

n (y)

=

• k,n

−k2uk,neikxTn(y) +

• k,n

uk,neikx

m

Rm,nTm(y) =

• k,n

−k2uk,neikxTn(y) +

• k,m
• n

Rm,nuk,n

• eikxTm(y)

=

• k,n

−k2uk,neikxTn(y) +

• k,m

(Ruk)meikxTm(y) =

• k,n

−k2uk,neikxTn(y) +

• k,n

(Ruk)neikxTn(y) =

• k,n
• −k2uk,n + (Ruk)n
• eikxTn(y)

where (Ruk)m ≡

m Rm,nuk,n

∆u(x, y) =

• k,n

 

k′,n′

∆k,n,k′,n′uk′,n′   eikxTn(y) where ∆k,n,k′,n′ = −k2δk,k′δn,n′ + Rn,n′δk,k′ Operators in x commute with operators in y. Odd and even Chebyshev polynomials are decoupled. Recall: second-derivative Chebyshev matrix R is upper triangular tridiagonal matrix B is such that BR is tridiagonal. B∆u(x, y) =

• k,n,n′
• −k2Bn,n′ + (BR)n,n′

uk′,n′eikxTn(y) Boundary conditions needed

    1 1 1 1 1 (BR − k2B)2,0 (BR − k2B)2,2 (BR − k2B)2,4 (BR − k2B)4,2 (BR − k2B)4,4 (BR − k2B)4,6 (BR − k2B)6,4 (BR − k2B)6,6 (BR − k2B)6,8    

Dirichlet BCs in x and y: Finite differences/Finite differences = Fill-in to maximal bandwidth: bandedness not preserved Bandwidth J = Nx, so operation count would be O(J2N) = O(N 2

xNxNy)

for LU decomposition and O(JN) = O(NxNxNy) for backsolve.

Alternative way of inverting: diagonalize in one or both directions. Recall that operators in x commute with operators in y so they are simultaneously diagonalizable ∆(n, m, n′, m′) = Dxx(n, n′)δ(m, m′) + Dyy(m, m′)δ(n, n′) =

• VxΛxV −1

x

• (n, n′)δ(m, m′) +
• VyΛyV −1

y

• (m, m′)δ(n, n′)

Diagonalization (cost N 3

x + N 3 y):

Dxx = VxΛxV −1

x

Dyy = VyΛyV −1

y

Vx, Vy are the Nx × Nx and Ny × Ny matrices of eigenvectors Dxx, Dyy are Nx × Nx and Ny × Ny matrices representing ∂xx, ∂yy Λx and Λy are diagonal matrices containing the eigenvalues Transform to x and y eigenspace O(N 2

xNy + N 2 yNx)

Invert Laplacian in eigenspace: O(NxNy) Inverse transform back from x and y eigenspace O(N 2

xNy + N 2 yNx)

Total: O(NxNy(Nx + Ny)) Storage: O(N 2

x + N 2 y)

Can use diag in x and banded LU in y with cost O(NxNy(Nx + 3)) Can also do in 3D, with timing O(NxNyNz(Nx + Ny + Nz))

Write f(x, y) as Nx × Ny matrix instead of vector of length NxNy ∂xxf(xi, yj) =

• k

Dxx(i, k)f(xk, yj) = (Dxxf)(xi, yj) ∂yyf(xi, yj) =

• k

Dyy(j, k)f(xi, yk) = (fDT

yy)(xi, yj)

∂xxf = Dxxf = VxΛxV −1

x

f ∂yyf = fDT

yy = f(VyΛyV −1 y

)T = f(V −1

y

)TΛyV T

y

∇2f = g VxΛxV −1

x

f + f(V −1

y

)TΛyV T

y

= g Λx

• V −1

x

f(V T

y )−1

+

• V −1

x

f(V −1

y

)T Λy =

• V −1

x

g(V T

y )−1

• V −1

x

f(V T

y )−1

(i, j) =

• V −1

x

g(V T

y )−1

(i, j) Λx(i) + Λy(j) f = Vx

• V −1

x

f(V T

y )−1

V T

y

Stencils for two-dimensional finite-difference Laplacian

Five-point stencil: d2u dx2(xj, yk) + d2u dy2(xj, yk) ≈ 1 h2 (u(xj+1, yk) − 2u(xj, yk) + u(xj−1, yk)) + 1 h2 (u(xj, yk+1) − 2u(xj, yk) + u(xj, yk−1)) Error is h2 24(∂4

x + ∂4 y)u + . . .

This error is not isotropic, unlike the Laplacian itself.

We show that the continuous Laplacian IS isotropic: Rotate (x, y) → (x′, y′) : x′ = αx+βy, y′ = −βx+αy, α2+β2 = 1 ∂ ∂x = ∂x′ ∂x ∂ ∂x′ + ∂y′ ∂x ∂ ∂y′ = α ∂ ∂x′ − β ∂ ∂y′ ∂2 ∂x2 =

• α ∂

∂x′ − β ∂ ∂y′ 2 = α2 ∂2 ∂x′2 − 2αβ ∂2 ∂x′∂y′ + β2 ∂2 ∂y′2 ∂ ∂y = ∂x′ ∂y ∂ ∂x′ + ∂y′ ∂y ∂ ∂y′ = β ∂ ∂x′ + α ∂ ∂y′ ∂2 ∂y2 =

• β ∂

∂x′ + α ∂ ∂y′ 2 = β2 ∂2 ∂x′2 + 2αβ ∂2 ∂x′∂y′ + α2 ∂2 ∂y′2 ∂2 ∂x2 + ∂2 ∂y2 = ∂2 ∂x′2 + ∂2 ∂y′2

Error of nine-point stencil is isotropic (to lowest order) d2u dx2(xj, yk) + d2u dy2(xj, yk) ≈ 1 6h2 [4 (u(xj+1, yk) + u(xj−1, yk) + u(xj, yk+1) + u(xj, yk−1)) + u(xj+1, yk+1) + u(xj−1, yk−1) + u(xj+1, yk−1) + u(xj−1, yk+1) −20u(xj, yk)] Error: h2 12

• ∂2

x + ∂2 y

2 = h2 12∆2

An important iterative method for solving the Poisson equation:

Multigrid

Like the FFT, the multigrid algorithm relies on coarsening the grid recursively, solving on coarse grids, then returning to fine grids. An additional reason to emphasize solutions to the Poisson equation:

Helmholtz problems

Parabolic problems (like heat equation, Navier-Stokes equations, ...) Implicit schemes lead to Helmholtz-type problems such as (I − ∆t∇2)f = g which can be solved by the same techniques as the Poisson problem.