Draft Computation of frequency responses of PDEs in Chebfun - PowerPoint PPT Presentation

Draft Computation of frequency responses of PDEs in Chebfun Mihailo Jovanovi´ c www.umn.edu/ ∼ mihailo joint work with: Binh K. Lieu 2012 SIAM Annual Meeting

Draft M. J OVANOVI ´ C , UMN 1 Example: heat equation • Distributed input and output fields ϕ t ( y, t ) = ϕ yy ( y, t ) + d ( y, t ) ϕ ( y, 0) = 0 ϕ ( ± 1 , t ) = 0 ⋆ Harmonic forcing steady-state response d ( y, t ) = d ( y, ω ) e j ωt ϕ ( y, t ) = ϕ ( y, ω ) e j ωt − − − − − − − − − − − − − − − − − − − → ⋆ Frequency response operator ϕ ( y, ω ) = [ T ( ω ) d ( · , ω ) ] ( y ) (j ωI − ∂ yy ) − 1 d ( · , ω ) � � = ( y ) � 1 = T ker ( y, η ; ω ) d ( η, ω ) d η − 1

Draft M. J OVANOVI ´ C , UMN 2 Two point boundary value realizations of T ( ω ) • Input-output differential equation � ϕ ′′ ( y, ω ) − j ω ϕ ( y, ω ) − d ( y, ω ) = T ( ω ) : ϕ ( ± 1 , ω ) = 0 • Spatial state-space realization � x ′ � 0 � � x 1 ( y, ω )  1 ( y, ω ) � 1 � � 0 � = + d ( y, ω )   − 1  x ′ 2 ( y, ω ) j ω 0 x 2 ( y, ω )   0 � � x 1 ( y, ω )   � 1  �  T ( ω ) : ϕ ( y, ω ) = x 2 ( y, ω )  � � x 1 ( − 1 , ω ) � � x 1 (1 , ω )   � � � �  1 0 0 0   0 = +   0 0 1 0 x 2 ( − 1 , ω ) x 2 (1 , ω ) 

Draft M. J OVANOVI ´ C , UMN 3 Frequency response operator • Evolution equation [ E φ t ( · , t )] ( y ) = [ F φ ( · , t )] ( y ) + [ G d ( · , t )] ( y ) , y ∈ [ a, b ] ϕ ( y, t ) = [ H φ ( · , t )] ( y ) , t ∈ [0 , ∞ ) ⋆ Spatial differential operators n ij f ij,k ( y ) d k � F = [ F ij ] = d y k k = 0 ⋆ Frequency response operator = H (j ω E − F ) − 1 G T

Draft M. J OVANOVI ´ C , UMN 4 Example: channel flow • Streamwise-constant fluctuations (1 /Re ) ∆ 2 � � � � � � � � ∆ 0 φ 1 t 0 φ 1 = F 21 0 I φ 2 t (1 /Re ) ∆ φ 2 ∆ = ∂ yy − k 2 Laplacian: z ∆ 2 = ∂ yyyy − 2 k 2 z ∂ yy + k 4 ”Square of Laplacian”: z F 21 = − j k z U ′ ( y ) Coupling:

Draft M. J OVANOVI ´ C , UMN 5 Singular value decomposition • Schmidt decomposition of a compact operator T ( ω ) : H in − → H out ∞ � ϕ ( y, ω ) = [ T ( ω ) d ( · , ω )] ( y ) = σ n ( ω ) u n ( y, ω ) � v n , d � n = 1 • Left and right singular functions [ T ( ω ) T ⋆ ( ω ) u n ( · , ω )] ( y ) = σ 2 n ( ω ) u n ( y, ω ) [ T ⋆ ( ω ) T ( ω ) v n ( · , ω )] ( y ) = σ 2 n ( ω ) v n ( y, ω ) { u n } orthonormal basis of H out { v n } orthonormal basis of H in

Draft M. J OVANOVI ´ C , UMN 6 • Right singular functions ⋆ identify input directions with simple responses σ 1 ( ω ) ≥ σ 2 ( ω ) ≥ · · · > 0 ∞ � ϕ ( ω ) = T ( ω ) d ( ω ) = σ n ( ω ) u n ( ω ) � v n ( ω ) , d ( ω ) � n = 1  � d ( ω ) = v m ( ω )  ϕ ( ω ) = σ m ( ω ) u m ( ω ) σ 1 ( ω ) : the largest amplification at any frequency

Draft M. J OVANOVI ´ C , UMN 7 Worst case amplification • H ∞ norm: an induced L 2 gain (of a system) � ∞ � ϕ ( t ) , ϕ ( t ) � d t sup output energy 0 G = = sup � ∞ input energy � d ( t ) , d ( t ) � d t 0 ω σ 2 1 ( T ( ω )) = sup 1 ( T ( ω )) σ 2

Draft M. J OVANOVI ´ C , UMN 8 Robustness interpretation ϕ d Nominal s ✲ ✲ Dynamics Γ ✛ modeling uncertainty (can be nonlinear or time-varying) small-gain theorem: stability for all Γ with γ 2 < 1 ⇔ σ 2 1 (Γ( ω )) ≤ γ 2 max G ω LARGE ⇒ small stability margins worst case amplification closely related to pseudospectra of linear operators

Draft M. J OVANOVI ´ C , UMN 9 Pseudo-spectral methods • MATLAB Differentiation Matrix Suite j ωI − D (2) � − 1 T ( ω ) = � ≈ N N • Advantages ⋆ superior accuracy compared to finite difference methods ⋆ ease-to-use M ATLAB codes • Disadvantages ⋆ ill-conditioning of high-order differentiation matrices ⋆ implementation of boundary conditions may be non-trivial Weideman & Reddy, ACM. TOMS. ’00

Draft M. J OVANOVI ´ C , UMN 10 Alternative method 1. Frequency response operator: two-point boundary value problem 2. Integral form of differential equations 3. State-of-the-art automatic spectral collocation techniques

Draft M. J OVANOVI ´ C , UMN 11 Advantages of Chebfun • Superior accuracy compared to currently available schemes • Avoids ill-conditioning of high-order differentiation matrices • Incorporates a wide range of boundary conditions • Easy-to-use MATLAB codes Lieu & Jovanovi´ c “Computation of frequency responses of linear time-invariant PDEs on a compact interval”, submitted to J. Comput. Phys., 2011 Also arXiv:1112.0579v1

Draft M. J OVANOVI ´ C , UMN 12 2D inertialess flow of viscoelastic fluids − ∇ p + (1 − β ) ∇ · τ + β ∇ 2 v + d 0 = 0 = ∇ · v τ t = ∇ v + ( ∇ v ) T − τ + We ( τ · ∇ ¯ v + ¯ τ · ∇ v τ · ∇ v ) T + ( τ · ∇ ¯ v ) T − v · ∇ ¯ + (¯ τ − ¯ v · ∇ τ ) We = polymer relaxation time characteristic flow time

Draft M. J OVANOVI ´ C , UMN 13 Largest singular value of T σ 1 ( T ) β = 0 . 5 k x = 1 ω = 0 N = 50 ( × ) 100 ( ◦ ) N = N = 200 (+) We

Draft M. J OVANOVI ´ C , UMN 14 σ 1 ( T ) β = 0 . 5 k x = 1 ω = 0 N = 50 ( × ) 100 ( ◦ ) N = N = 200 (+) We

Draft M. J OVANOVI ´ C , UMN 15 Input-output differential equations • Frequency response operator  [ A 0 φ ] ( y ) = [ B 0 d ] ( y )    T ( ω ) : ϕ ( y ) = [ C 0 φ ] ( y )  0 = N 0 φ ( y )   • Adjoint of the frequency response operator  [ A ⋆ 0 ψ ] ( y ) = [ C ⋆ 0 f ] ( y )    T ⋆ ( ω ) : g ( y ) = [ B ⋆ 0 ψ ] ( y )  0 = N ⋆ 0 ψ ( y )  

Draft M. J OVANOVI ´ C , UMN 16 Composition of T with T ⋆ • Cascade connection  [ A ξ ] ( y ) = [ B f ] ( y )   T T ⋆ : ϕ ( y ) = [ C ξ ] ( y )  0 = N ξ ( y )  ⋆ Do e-value decomposition of T T ⋆ and T ⋆ T in Chebfun [ T ( ω ) T ⋆ ( ω ) u n ( · , ω )] ( y ) = σ 2 n ( ω ) u n ( y, ω ) [ T ⋆ ( ω ) T ( ω ) v n ( · , ω )] ( y ) = σ 2 n ( ω ) v n ( y, ω )

Draft M. J OVANOVI ´ C , UMN 17 Example: 1D heat equation y ∈ [ − 1 , 1] φ t ( y, t ) = φ yy ( y, t ) + d ( y, t ) , φ ( ± 1 , t ) = 0 • Frequency response and adjoint operators  φ ′′ ( y ) − j ω φ ( y ) = − d ( y )  �� 1 � 0 � 0  T ( ω ) : � � � � E − 1 + E 1 φ ( y ) =  0 1 0   ψ ′′ ( y ) + j ω ψ ( y ) = f ( y )     g ( y ) = − ψ ( y )  T ⋆ ( ω ) : �� 1 � 0 � 0 � � � �  E − 1 + E 1 ψ ( y ) =    0 1 0 

Draft M. J OVANOVI ´ C , UMN 18 Integral form of a differential equation Driscoll, J. Comput. Phys., 2010 • 1D diffusion equation: differential form � D (2) − j ωI � φ ( y ) = − d ( y ) �� � � � � � � 1 0 0 E − 1 + E 1 φ ( y ) = 0 1 0 D (2) φ � � Auxiliary variable: ν ( y ) = ( y ) Integrate twice � y � J (1) ν � φ ′ ( y ) = ν ( η 1 ) d η 1 + k 1 = ( y ) + k 1 − 1 � y �� η 2 � 1 � ( y + 1) � � � k 2 φ ( y ) = ν ( η 1 ) d η 1 d η 2 + k 1 − 1 − 1 J (2) ν ( y ) + K (2) k � � =

Draft M. J OVANOVI ´ C , UMN 19 • 1D diffusion equation: integral form � I − j ωJ (2) � ν ( y ) − j ω K (2) k = − d ( y ) � � � � �� � � � � 1 0 k 2 1 0 J (2) ν ( y ) = − E − 1 + E 1 1 2 k 1 0 1 Eliminate k from the equations to obtain � I − j ωJ (2) + 1 � 2 j ω ( y + 1) E 1 J (2) ν ( y ) = − d ( y ) ☞ More suitable for numerical computations than differential form integral operators and point evaluation functionals are well-conditioned

Draft M. J OVANOVI ´ C , UMN 20 Online resources • Chebfun http://www2.maths.ox.ac.uk/chebfun/ Google: “chebfun” • Computing frequency responses of PDEs http://www.umn.edu/ ∼ mihailo/software/chebfun-svd/ Google: “frequency responses pde”

Draft M. J OVANOVI ´ C , UMN 21 Summary • method for computing frequency responses of PDEs • easy-to-use mini-toolbox in M ATLAB ⋆ enabling tool: Chebfun • two major advantages over currently available schemes ⋆ avoids ill-conditioning of high-order differentiation matrices ⋆ easy implementation of boundary conditions

Draft M. J OVANOVI ´ C , UMN 22 Acknowledgments T EAM : Binh Lieu S UPPORT : NSF CAREER Award CMMI-06-44793 S OFTWARE : http://www.umn.edu/ ∼ mihailo/software/chebfun-svd/