u + u = 0 , 2 u ( ) = u ( ) = 1 The basic idea: represent the - PowerPoint PPT Presentation

Lesson 11 S PECTRAL METHODS

• In this lecture we consider the problem of computing solutions to linear ODEs: a ( N ) ( θ ) u ( N ) + · · · + a (0) ( θ ) u = f ( θ ) on the periodic interval a ( N ) ( z ) u ( N ) + · · · + a (0) ( z ) u = f ( z ) inside the unit disk a ( N ) ( x ) u ( N ) + · · · + a (0) ( x ) u = f ( x ) on the unit interval • We will impose boundary conditions so that the solution exists and is unique • For example, on the periodic interval we might solve: u �� + u ��� θ = 0 , 2 u ( − π ) = u ( π ) = 1 • The basic idea: represent the infinite-dimensional differential operator with a finite dimensional matrix

F OURIER SPECTRAL METHODS

• We start with a simple equation: u � + u = f ( θ ) , u ( − π ) = u ( π ) • Idea: represent u and f by their Fourier series: ∞ ∞ u k (i k + 1)e i k θ = ˆ � � f k e i k θ ˆ k =0 k =0 • The solution is obvious. But let’s anyways represent this in matrix form: .   .     ... . . . .   ˆ     f − 2 u − 2 − 2i + 1 ˆ         ˆ     f − 1 u − 1 − i + 1 ˆ             ˆ u 0 1 ˆ = f 0             ˆ u 1 i + 1 ˆ   f 1           u 2 2i + 1 ˆ ˆ   f 2           . ...  .  . . . .

• Numerically, we truncate the operator and use the DFT to compute the coefficients: � � α � + 1 � � ... � � u m = F f m � β � + 1 where � � f ( θ 1 ) . � � . � , and f m = f ( θ m ) = � . f ( θ m ) � � ˆ u α . � � . � u m ≈ � � . ˆ u β • In this simple case, convergence follows from convergence of f m

• An alternative is to represent the operator as acting on function values:   α � + 1 ... F − 1  F u m = f m    β � + 1 where u m ≈ u ( θ m ) • Which method is better? Let’s compare

Matrix density coefficient space value space O( m ) operations O( m 3 ) operations

� • We introduce some notation • The differentiation operator is the bi-infinite operator � � ... � � � � � � − � � � D = 0 � � � � � � ... • The truncation operator maps bi-infinite vectors to finite vectors: � � u α . � � . P m � u = � � . u β • The truncated differentiation matrix is: � � � α � � ... D m = P m DP � m = � � . � β

• We want to now incorporate variable coefficients, i.e., equations of the form u � + a ( θ ) u = f ( θ ) • In value space this is: � � � � a ( θ 1 ) � � � � ... � F � 1 D m F + � u m = f m � � a ( θ m ) where D m is the discretization of the differentiation operator: � � � α � � ... D m = � . � � β • In coefficient space we have: � � � � a ( θ 1 ) � � � � ... � F � 1 � D m + F � � u m = F f m � a ( θ m )

u � ( θ ) + u ( θ ) ��� θ = 1 u ( − π ) = u ( π ) 1 0.001 10 - 6 Error 10 - 9 value space 10 - 12 10 - 15 coefficient space 50 100 150 200 250 300 350 Number of unknowns

Matrix density coefficient space value space O( C m ) operations O( m 3 ) operations

M ULTIPLICATION OF F OURIER SERIES

• We want to simplify � � a ( θ 1 ) ... � F − 1 , F � � � a ( θ m ) i.e., we want to investigate multiplication in coefficient space • Let's consider a ( θ ) = � − � θ : we have ∞ ∞ f k � � k θ = � ˆ � ˆ � − � θ f k +1 � � k θ k = −∞ k = −∞ • Thus in operator space we have the shift operator: . . ... . . . . so that ... . . . . . .

• We want to simplify � � a ( θ 1 ) ... � F − 1 , F � � � a ( θ m ) i.e., we want to investigate multiplication in coefficient space • Let's consider a ( θ ) = � − � θ : we have ∞ ∞ f k � � k θ = � ˆ � ˆ � − � θ f k +1 � � k θ k = −∞ k = −∞ • Thus in operator space we have the shift operator: . . � � � � ... . . � � . . � � � � ˆ ˆ � � f − 1 f 0 1 � � � � � � � � � � ˆ ˆ � � , so that S 0 1 S = = f 0 f 1 � � � � � � � � � � ... � � ˆ ˆ � � � � f 1 f 2 � � � � � � � � . . � � � � . . . .

� � � • More generally, if a ( θ ) = � � j θ : we have ∞ ∞ � � f k � � k θ = ˆ ˆ � − � j θ f k + j � � k θ k = −∞ k = −∞ times � , i.e., • This is equivalent to ... ...

• More generally, if a ( θ ) = � � j θ : we have ∞ ∞ � � f k � � k θ = ˆ ˆ � − � j θ f k + j � � k θ k = −∞ k = −∞ j times � �� � � − � θ · · · � − � θ , i.e., • This is equivalent to � � ... � � 1 � � � � S j = 0 0 0 1 · · · � � � � ... � � j times � �

• For general, smooth and periodic a we can expand in Fourier series: ∞ � a k � � k θ a ( θ ) = ˆ k = −∞ • Each term is multiplication by � � j θ • Thus, in operator space we have the multiplication operator ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... • A bi-infinite matrix with constant diagonals is called a Laurent operator

• For general, smooth and periodic a we can expand in Fourier series: ∞ � a k � � k θ a ( θ ) = ˆ k = −∞ • Each term is multiplication by � � j θ • Thus, in operator space we have the multiplication operator ... ... ... ... ... � � ... ... � � � � ˆ ˆ ˆ a 0 a − 1 a − 2 � � ∞ ... ... � � a k S − k = � L [ a ] = ˆ � � ˆ ˆ ˆ a 1 a 0 a − 1 � � ... ... � � k = −∞ � � ˆ ˆ ˆ a 2 a 1 a 0 � � � ... ... ... ... ... � • A bi-infinite matrix with constant diagonals is called a Laurent operator

• In practice, we can approximate a : β � k � � k θ a m a ( θ ) ≈ ˆ k = α • Thus, in operator space the multiplication operator is banded • For example, if we have

• In practice, we can approximate a : β � k � � k θ a m a ( θ ) ≈ ˆ k = α • Thus, in operator space the multiplication operator is banded • For example, if α = β = 1 we have ... ... ... ... ... � � ... ... � � � � ... ... � � ... � � ˆ ˆ ˆ a 0 a − 1 a − 2 � a m a m � ˆ ˆ � � � � 0 − 1 � � ... ... � a m a m a m � ˆ ˆ ˆ L [ a ] = � � ≈ ˆ ˆ ˆ a 1 a 0 a − 1 � � 1 0 − 1 � � � � ... � � ... ... � � a m a m ˆ ˆ � � ˆ ˆ ˆ � � a 2 a 1 a 0 1 0 � � � � ... ... � � ... ... ... ... ...

• Finally, we truncate the multiplication operator. Thus the operator u � + a ( θ ) u = f becomes ( D m + L m [ a ]) � u m = F f m where � � ˆ ˆ a 0 a � 1 � � ˆ ˆ ˆ a 1 a 0 a � 1 � � � � ... ... ... L m [ a ] = P m L [ a ] P � m = � � � � � � ˆ ˆ ˆ a 1 a 0 a � 1 ˆ ˆ a 1 a 0

H IGHER ORDER EQUATIONS

• The operator a ( N ) ( θ ) u ( N ) + · · · + a (0) u in coefficient space is L [ a ( N ) ] D N + · · · L [ a (0) ] • We can truncate it by L m [ a ( N ) ] D N m + · · · L m [ a (0) ] • However, we need to impose boundary conditions so that the operator is uniquely invertible • In finite dimensions, non-unique inverse causes bad conditioning

• We denote the boundary conditions as an d × ∞ operator: B � u • An example: suppose we want to impose Neumann conditions at the point χ : u ( χ ) = u � ( χ ) = 0 � Then we have the boundary condition operator � � � � � 2 χ � 2 � χ � � � χ 1 · · · · · · � � χ B = − 2 �� � � 2 χ 2 �� 2 � χ − �� � � χ 0 · · · · · · �� � χ

• The differential equation a ( N ) ( θ ) u ( N ) + · · · + a (0) ( θ ) u = f ( θ ) and B u = c in coefficient space is a pair of equations: and B � u = c � � L [ a ( N ) ] D N + · · · + L [ a (0) ] u = � � f • We need a square matrix, so we now truncate it as � � � � BP � c � � m u m = � L [ a ( N ) ] D N + · · · + L [ a (0) ] P � F f m � d P m � d m

T AYLOR SERIES SPECTRAL METHODS

• For Taylor series, the vector of coefficients becomes only infinite in one direction � � ˆ u 0 � � ˆ u 1 � u = � � . . . • Differentiation becomes � � 0 1 � � 2 � � � � D = 3 � � ... • Truncation is now defined by � � ˆ u 0 . � � . P m � u = � � . ˆ u m − 1