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 N − 1 � � u ( x ℓ ) e − i ( k + N )2 πℓ/N = u ( x ℓ ) e − ik 2 πℓ/N e − iN 2 πℓ/N u k + N = ˆ ℓ =0 � u ( x ℓ ) e − ik 2 πℓ/N = ˆ = u k u k is N -periodic ˆ N/ 2 − 1 N/ 2 − 1 � � u ( x 2 ℓ ) e − i 2 π (2 ℓ ) k/N + u ( x 2 ℓ +1 ) e − i 2 π (2 ℓ +1) k/N u k = ˆ ℓ =0 ℓ =0 N/ 2 − 1 N/ 2 − 1 � � u ( x 2 ℓ ) e − i 2 πℓk/ ( N/ 2) + e − i 2 πk/N u ( x 2 ℓ +1 ) e − i 2 πℓk/ ( N/ 2) = ℓ =0 ℓ =0 N/ 2 − 1 N/ 2 − 1 � � v ( x ℓ ) e − i 2 πℓk/ ( N/ 2) + e − i 2 πk/N w ( x ℓ ) e − i 2 πℓk/ ( N/ 2) = ℓ =0 ℓ =0 + e − i 2 πk/N = v k ˆ w k ˆ − e − i 2 πk/N u k + N/ 2 = ˆ v k ˆ w k ˆ 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 + N = N 2 N 2 − → 2 + N 2 2 � � � N � 2 + N = N 2 N + N 4 − → 2 2 + 2 N 4 2 4 � � � � � N � 2 + N = N 2 N + N + N 8 − → 2 2 2 + 3 N 8 4 2 8 N = 2 p : → N 2 N 2 p − 2 p + pN = N + (log 2 N ) N = O ( N log 2 N )

Fourier transform in one dimension: N − 1 � u ( x ℓ ) e − i 2 πkℓ/N u k = � ℓ =0 One multiplication for each k, ℓ , so O ( N 2 x ) operations. Fourier transform in two dimensions: � N x − 1 � N y − 1 � � � u ( x ℓ , y n ) e − i 2 πkℓ/N x e − i 2 πmn/N y � u k,m = n =0 ℓ =0 � �� � u k ( y n ) � N x N 2 N 2 x N y { � y { � � { u ( x ℓ , y n ) } − − − → u k ( y n ) } − − − → u k,m } Total: N x N y ( N x + N y ) ≪ N 2 x N 2 y With FFT: N x N y (log N x + log N y ) = N x N y log( N x N y )

Even without FFT, multidimensional Fourier transform would be fast be- cause the different dimensions are decoupled: N x N y N z ( N x + N y + N z ) 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 � � � � � � f ∗ � f ℓ � fg ( k ) = g ( k ) = g k − ℓ ℓ One multiplication for each k, ℓ , so O ( N 2 x ) operations. Convolution in two dimensions: � � � � � � � f ∗ � f ℓ,n � fg ( k, m ) = g ( k, m ) = g k − ℓ,m − n n ℓ One multiplication for each k, ℓ, m, n , so O ( N 2 x N 2 y ) operations.

Multidimensional Fourier transform of real data u ∗ For u ( x ) real, ˆ u k is conjugate symmetric: ˆ u − k = ˆ k so that � u k e ikx � ∗ = 2 R e � u k e ikx � u k e ikx + ˆ u − k e − ikx = ˆ u k e ikx + u ( x ) ∼ ˆ ˆ ˆ u ∗ For u ( x, y ) real, ˆ u − k, − m = ˆ k,m so need half of ( k, m ) plane. or u ∗ For u ( x, y, z ) real, ˆ u − k, − m, − n = ˆ 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 k,m e ikx e imy g k,m e ikx e imy u ( x, y ) = ˆ g ( x, y ) = ˆ k,m k,m � � ( − k 2 − m 2 )ˆ u k,m e ikx e imy = g k,m e ikx e imy ∆ u = ˆ k,m k,m − ˆ g k,m u k,m = ˆ k 2 + m 2 Have used interval [0 , 2 π ) for simplicity. More generally, domain is [0 , L x ) × [0 , L y ) and basis functions are e ikx 2 π/L x e imy 2 π/L y

Periodic BCs in x , Dirichlet BCs in y : Fourier-Finite Differences � � u k ( y ) e ikx g k ( y ) e ikx u ( x, y ) = ˆ g ( x, y ) = ˆ k k � � � u k ( y ) + ˆ u k ( y + ∆ y ) − 2ˆ u k ( y ) + ˆ u k ( y − ∆ y ) − k 2 ˆ e ikx ∆ u = (∆ y ) 2 k � − 2 − ( k ∆ y ) 2 � � u k ( y + ∆ y ) + ˆ u k ( y ) + ˆ ˆ u k ( y − ∆ y ) e ikx = (∆ y ) 2 k For each Fourier component k     − 2 − ( k ∆ y ) 2 1 u k ( y 1 ) ˆ     − 2 − ( k ∆ y ) 2 1 1 u k ( y 2 ) ˆ     1     − 2 − ( k ∆ y ) 2 1 1 u k ( y 3 ) ˆ         (∆ y ) 2 . ... .     . − 2 − ( k ∆ y ) 2 ) 1 u k ( y N ) ˆ

Boundary conditions needed for each Fourier component k , e.g.       1 0 0 . . . 0 u k ( y 1 ) ˆ γ −       − 2 − ( k ∆ y ) 2 1 1 u k ( y 2 ) ˆ g k ( y 2 ) ˆ       1       − 2 − ( k ∆ y ) 2 1 1 u k ( y 3 ) ˆ ˆ g k ( y 3 )     =         (∆ y ) 2 . . ... . .       . . 0 0 0 . . . 1 u k ( y N ) ˆ γ + Elliptic equations need boundary conditions (Dirichlet, Neumann or peri- odic) along all boundaries of domain

Periodic BCs in x , Dirichlet BCs in y : Fourier-Chebyshev � � u k,n e ikx T n ( y ) g k,n e ikx T n ( y ) u ( x, y ) = g ( x, y ) = k,n k,n � � − k 2 u k,n e ikx T n ( y ) + u k,n e ikx T ′′ ∆ u = n ( y ) k,n k,n � � u k,n e ikx � − k 2 u k,n e ikx T n ( y ) + = R m,n T m ( y ) m k,n k,n �� � � � − k 2 u k,n e ikx T n ( y ) + e ikx T m ( y ) = R m,n u k,n n k,n k,m � � − k 2 u k,n e ikx T n ( y ) + ( Ru k ) m e ikx T m ( y ) = k,n k,m � � − k 2 u k,n e ikx T n ( y ) + ( Ru k ) n e ikx T n ( y ) = k,n k,n � � � − k 2 u k,n + ( Ru k ) n e ikx T n ( y ) = k,n where ( Ru k ) m ≡ � m R m,n u k,n

  � �  e ikx T n ( y ) ∆ u ( x, y ) = ∆ k,n,k ′ ,n ′ u k ′ ,n ′ k,n k ′ ,n ′ ∆ k,n,k ′ ,n ′ = − k 2 δ k,k ′ δ n,n ′ + R n,n ′ δ k,k ′ where 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. � � − k 2 B n,n ′ + ( BR ) n,n ′ � u k ′ ,n ′ e ikx T n ( y ) B ∆ u ( x, y ) = k,n,n ′ Boundary conditions needed   1 1 1 1 1   ( BR − k 2 B ) 2 , 0 ( BR − k 2 B ) 2 , 2 ( BR − k 2 B ) 2 , 4     ( BR − k 2 B ) 4 , 2 ( BR − k 2 B ) 4 , 4 ( BR − k 2 B ) 4 , 6 ( BR − k 2 B ) 6 , 4 ( BR − k 2 B ) 6 , 6 ( BR − k 2 B ) 6 , 8

Dirichlet BCs in x and y : Finite differences/Finite differences = Fill-in to maximal bandwidth: bandedness not preserved Bandwidth J = N x , so operation count would be O ( J 2 N ) = O ( N 2 x N x N y ) for LU decomposition and O ( JN ) = O ( N x N x N y ) 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 ′ ) = D xx ( n, n ′ ) δ ( m, m ′ ) + D yy ( m, m ′ ) δ ( n, n ′ ) � � � � V x Λ x V − 1 ( n, n ′ ) δ ( m, m ′ ) + V y Λ y V − 1 ( m, m ′ ) δ ( n, n ′ ) = x y Diagonalization (cost N 3 x + N 3 y ): D xx = V x Λ x V − 1 D yy = V y Λ y V − 1 x y V x , V y are the N x × N x and N y × N y matrices of eigenvectors D xx , D yy are N x × N x and N y × N y matrices representing ∂ xx , ∂ yy Λ x and Λ y are diagonal matrices containing the eigenvalues Transform to x and y eigenspace O ( N 2 x N y + N 2 y N x ) Invert Laplacian in eigenspace: O ( N x N y ) Inverse transform back from x and y eigenspace O ( N 2 x N y + N 2 y N x ) Storage: O ( N 2 x + N 2 Total: O ( N x N y ( N x + N y )) y ) Can use diag in x and banded LU in y with cost O ( N x N y ( N x + 3)) Can also do in 3D, with timing O ( N x N y N z ( N x + N y + N z ))

Write f ( x, y ) as N x × N y matrix instead of vector of length N x N y � ∂ xx f ( x i , y j ) = D xx ( i, k ) f ( x k , y j ) = ( D xx f )( x i , y j ) k � D yy ( j, k ) f ( x i , y k ) = ( fD T ∂ yy f ( x i , y j ) = yy )( x i , y j ) k = D xx f = V x Λ x V − 1 ∂ xx f f x ) T = f ( V − 1 = fD T yy = f ( V y Λ y V − 1 ) T Λ y V T ∂ yy f y y y ∇ 2 f = g V x Λ x V − 1 f + f ( V − 1 ) T Λ y V T = g x y y � y ) − 1 � � ) T � � y ) − 1 � V − 1 f ( V T V − 1 f ( V − 1 V − 1 g ( V T Λ x + Λ y = x x y x � y ) − 1 � V − 1 g ( V T � y ) − 1 � ( i, j ) x V − 1 f ( V T ( i, j ) = x Λ x ( i ) + Λ y ( j ) � y ) − 1 � V − 1 f ( V T V T f = V x x y

Stencils for two-dimensional finite-difference Laplacian Five-point stencil: d 2 u dx 2 ( x j , y k ) + d 2 u dy 2 ( x j , y k ) ≈ 1 h 2 ( u ( x j +1 , y k ) − 2 u ( x j , y k ) + u ( x j − 1 , y k )) + 1 h 2 ( u ( x j , y k +1 ) − 2 u ( x j , y k ) + u ( x j , y k − 1 )) h 2 24( ∂ 4 x + ∂ 4 y ) u + . . . Error is This error is not isotropic, unlike the Laplacian itself.