Draft EE 8235: Lectures 17 & 18 1 Lectures 17 & 18: - PowerPoint PPT Presentation

Draft EE 8235: Lectures 17 & 18 1 Lectures 17 & 18: Numerical methods • Spectral (Galerkin) method ⋆ Basis function expansion ⋆ Compute inner products to determine equation for spectral coefficients • Pseudo-spectral method ⋆ Satisfy equation at the set of ”collocation” points ⋆ Connection to polynomial interpolation • Chebyshev polynomials ⋆ Why they should be used ⋆ Basic properties

Draft EE 8235: Lectures 17 & 18 2 Online resources • Freely available books/papers ⋆ Jonh P . Boyd Chebyshev and Fourier Spectral Methods ⋆ Lloyd N. Trefethen Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations ⋆ Weideman and Reddy A Matlab Differentiation Matrix Suite • Publicly available software ⋆ A Matlab Differentiation Matrix Suite http://dip.sun.ac.za/ ∼ weideman/research/differ.html ⋆ Chebfun http://www2.maths.ox.ac.uk/chebfun/

Draft EE 8235: Lectures 17 & 18 3 Diffusion equation on L 2 [ − 1 , 1] ψ t ( x, t ) = ψ xx ( x, t ) ψ ( x, 0) = ψ 0 ( x ) ψ ( ± 1 , t ) = 0 Basis function expansion ∞ � ψ ( x, t ) = α n ( t ) φ n ( x ) n = 1 α n ( t ) − (unknown) spectral coefficients φ n ( x ) − (known) basis functions

Draft EE 8235: Lectures 17 & 18 4 Galerkin method • Approximate solution by   α 1 ( t ) N � φ 1 ( x ) � . φ N ( x ) �   . ψ ( x, t ) ≈ α n ( t ) φ n ( x ) = · · · .   n = 1 α N ( t ) substitute into the equation and take an inner product with { φ m }         � φ 1 , φ ′′ � φ 1 , φ ′′ � φ 1 , φ 1 � · · · � φ 1 , φ N � α 1 ( t ) ˙ 1 � · · · N � α 1 ( t ) . . . . . .         . . . . . .  = . . . . . .        � φ N , φ ′′ � φ N , φ ′′ � φ N , φ 1 � · · · � φ N , φ N � α N ( t ) ˙ 1 � · · · N � α N ( t ) • Done if basis functions satisfy BCs Otherwise, need additional conditions on spectral coefficients   α 1 ( t ) � � � � 0 φ 1 ( − 1) · · · φ N ( − 1) .   . = .   0 φ 1 (+1) · · · φ N (+1) α N ( t )

Draft EE 8235: Lectures 17 & 18 5 Pros and cons • Advantage: superior convergence (if basis functions selected properly) • Problem: requires integration ⋆ Cumbersome in spatially-varying and nonlinear systems Example: Orr-Sommerfeld equation in fluid mechanics � � j k x ( U ′′ ( y ) − U ( y ) ∆) + 1 R ∆ 2 ∆ ψ t = ψ

Draft EE 8235: Lectures 17 & 18 6 Polynomial interpolation • Approximate f ( x ) by a polynomial that matches f ( x ) at interpolation points p N − 1 ( x i ) = f ( x i ) , i = { 1 , . . . , N } • Examples: N = 2 ⇒ Linear Interpolation N = 3 ⇒ Quadratic Interpolation x 2 x 0 x 1 x 1 x 0 ( x − x 1 )( x − x 2 ) f ( x ) ≈ ( x 0 − x 1 )( x 0 − x 2 ) f ( x 0 ) + f ( x ) ≈ x − x 1 f ( x 0 ) + x − x 0 ( x − x 0 )( x − x 2 ) f ( x 1 ) ( x 1 − x 0 )( x 1 − x 2 ) f ( x 1 ) + x 0 − x 1 x 1 − x 0 ( x − x 0 )( x − x 1 ) ( x 2 − x 0 )( x 2 − x 1 ) f ( x 2 )

Draft EE 8235: Lectures 17 & 18 7 Lagrange interpolation formula N � p N ( x ) = f ( x i ) C i ( x ) i = 0 N x − x j � C i ( x ) = x i − x j j = 0 , j � = i • Cardinal functions C i ( x j ) = δ ij ⋆ Not efficient for computations ⋆ Suitable for theoretical arguments • Runge Phenomenon 1 f ( x ) = 1 + x 2 , x ∈ [ − 5 , 5] ⋆ Evenly spaced points ⇒ convergence for | x | ≤ 3 . 63 Interactive Demo

Draft EE 8235: Lectures 17 & 18 8 Choice of grid points • Cauchy interpolation error theorem � − has N + 1 derivatives N f ⇒ f ( x ) − p N ( x ) = f ( N +1) ( ξ ) � ( x − x i ) ( N + 1)! p N − interpolant of degree N i = 0 • Chebyshev minimal amplitude theorem ⋆ Among all polynomials q N ( x ) of degree N , with leading coefficient 1 , T N ( x ) = N th Chebyshev polynomial 2 N − 1 2 N − 1 has the smallest L ∞ [ − 1 , 1] norm � � T N ( x ) 1 � � sup | q N ( x ) | ≥ sup � = for all q N ( x ) 2 N − 1 , � � 2 N − 1 � x ∈ [ − 1 , 1] x ∈ [ − 1 , 1]

Draft EE 8235: Lectures 17 & 18 9 Optimal interpolation points • Select polynomial part of f ( x ) − p N ( x ) as N ( x − x i ) = T N +1 ( x ) � 2 N i = 0 • Optimal interpolation points: roots of T N +1 ( x ) � (2 i − 1) π � x i = cos , i = { 1 , . . . , N + 1 } 2 ( N + 1)

Draft EE 8235: Lectures 17 & 18 10 Chebyshev polynomials • Solutions to Sturm-Liouville Problem � 1 − x 2 � n ( x ) + n 2 T n ( x ) = 0 , x ∈ [ − 1 , 1] , n = 0 , 1 , . . . T ′′ n ( x ) − x T ′ • Three-term recurrence { T 0 = 1; T 1 ( x ) = x ; T n +1 ( x ) = 2 x T n ( x ) − T n − 1 ( x ) , n ≥ 1 } • Alternative definition T n (cos ( t )) = cos ( n t ) ⇒ | T n ( x ) | ≤ 1 , for all x ∈ [ − 1 , 1] , n = 0 , 1 , . . .

Draft EE 8235: Lectures 17 & 18 11 • Inner product  0 m � = n   � 1   T m ( x ) T n ( x ) m = n = 0 π √ � T m , T n � w = d x = 1 − x 2  − 1 π   m = n � = 0  2 • Collocation points � (2 i − 1) π � Gauss-Chebyshev: x i = cos i = { 1 , . . . , N } , 2 N � � π i Gauss-Lobatto: x i = cos i = { 0 , . . . , N − 1 } , N − 1 • Integration � x T n +1 ( x ) T n − 1 ( x ) T n ( ξ ) d ξ = 2 ( n + 1) + 2 ( n − 1) , n ≥ 2 − 1

Draft EE 8235: Lectures 17 & 18 12 Gaussian integration • Approximate f ( x ) by a polynomial that matches f ( x ) at interpolation points p N ( x i ) = f ( x i ) , i = { 0 , . . . , N } N � f ( x ) ≈ p N ( x ) = f ( x i ) C i ( x ) i = 0 • Evaluate integral of f ( x ) by integrating p N ( x ) � b N � f ( x ) d x ≈ w i f ( x i ) a i = 0 Quadrature weights: � b w i = C i ( x ) d x a • Gaussian integration: exact if integrand is a polynomial of degree N

Draft EE 8235: Lectures 17 & 18 13 • Can be made exact for polynomials of degree 2 N + 1 by optimal selection of ⋆ interpolation points { x i } ⋆ weights { w i } • Gauss-Jacobi integration ⋆ orthogonal polynomials w.r.t. the inner product with weight function ρ ( x ) ⋆ interpolation points: zeros of p N +1 ( x ) ⋆ quadrature formula: exact for polynomials of degree 2 N + 1 or smaller � b N � f ( x ) ρ ( x ) d x = w i f ( x i ) a i = 0 • Good candidates for quadrature points: � π i � Gauss-Lobatto: x i = cos i = { 0 , . . . , N } , N

Draft EE 8235: Lectures 17 & 18 14 Interpolation by quadrature • Orthogonality w.r.t. discrete inner product N � � φ i , φ j � = δ ij ⇒ � φ i , φ j � G = w m φ i ( x m ) φ j ( x m ) = δ ij m = 0 • Basis function expansion ∞ N � � f ( x ) = α n φ n ( x ) = α n φ n ( x ) + E N ( x ) n = 0 n = 0 • Discrete vs. exact spectral coefficients = � φ m , f � G α m,G � � N � = α n φ n + E N φ m , n = 0 G N � = α n � φ m , φ n � G + � φ m , E N � G n = 0 = α m + � φ m , E N � G

Draft EE 8235: Lectures 17 & 18 15 Error bound for Chebyshev interpolation • Error between Galerkin and Pseudo-spectral twice the sum of absolute values of neglected spectral coefficients ∞ � ⋆ f ( x ) = α n T n ( x ) n = 0 ⋆ p N ( x ) – polynomial that interpolates f ( x ) at Gauss-Lobatto points ∞ � | f ( x ) − p N ( x ) | ≤ 2 | α n | , for all N and all x ∈ [ − 1 , 1] n = N +1

Draft EE 8235: Lectures 17 & 18 16 Back to cardinal functions • Lagrange interpolation N � p N ( x ) = f ( x i ) C i ( x ) i = 0 N x − x j � C i ( x ) = x i − x j j = 0 , j � = i Cardinal functions C i ( x j ) = δ ij • Sinc functions � ( x − kh ) π � sin � x − kh � h C k ( x ; h ) = = sinc ( x − kh ) π h h { x j = j h ; j ∈ Z } ⇒ C k ( x j ; h ) = δ jk Approximate f by ∞ � f ( x ) = f ( x j ) C j ( x ; h ) j = −∞

Draft EE 8235: Lectures 17 & 18 17 Cardinal functions for Chebyshev polynomials • Gauss-Chebyshev points: zeros of T N +1 ( x ) ⋆ Taylor series expansion around x j N +1 ( x j ) ( x − x j ) + 1 N +1 ( x j ) ( x − x j ) 2 + O � | x − x j | 3 � + T ′ 2 T ′′ T N +1 ( x ) = T N +1 ( x j ) ) � �� � 0 Cardinal functions N +1 ( x j ) ( x − x j ) = 1 + T ′′ N +1 ( x j ) ( x − x j ) T N +1 ( x ) � | x − x j | 2 � C j ( x ) = + O ) T ′ 2 T ′ N +1 ( x j ) • Gauss-Lobatto points: zeros of (1 − x 2 ) T ′ N ( x ) � 1 − x 2 � T ′ N ( x ) C j ( x ) = Cardinal functions: N ( x )) ′ � ((1 − x 2 ) T ′ x = x j ( x − x j ) �

Draft EE 8235: Lectures 17 & 18 18 Matlab Differentiation Matrix Suite: A Demo %% number of grid points without boundaries (no \pm 1) N = 50 %% 1st & 2nd order differentiation matrices [yT,DM] = chebdif(N+2,2); y = yT(2:end-1); %% 1st & 2nd derivatives wrt y on a total grid (no BCs) DT1 = DM(:,:,1); DT2 = DM(:,:,2); %% implement homogeneous Dirichlet BCs %% ammounts to deleting 1st rows and columns of DT1 & DT2 D1 = DT1(2:N+1,2:N+1); D2 = DT2(2:N+1,2:N+1); %% 4th derivative with Dirichlet & Neumann BCs at both ends %% D4 - obtained on a grid without \pm 1 [y1,D4] = cheb4c(N+2); %% e-value decomposition of D2 with Dirichlet BCs [Vh,Dh] = eig(D2); % compare with analytical results