Laurette TUCKERMAN laurette@pmmh.espci.fr Numerical Methods for Differential Equations in Physics
Curvilinear Coordinates Cylindrical coordinates ( r, θ, z ) (sometimes called ( ρ, θ, z ) ) r 2 = x 2 + y 2 x = r cos θ y = r sin θ θ = atan2 ( y, x ) z = z z = z Spherical coordinates ( r, θ, φ ) (sometimes called ( r, φ, θ ) ) r 2 = x 2 + y 2 + z 2 x = r sin θ cos φ � x 2 + y 2 , z ) y = r sin θ sin φ θ = atan2 ( z = r cos θ φ = atan2 ( y, x )
� � − π 2 , π Range of atan is 2 Range of atan2 is ( − π, π ] Differential operators in cylindrical coordinates contain 1 r ∇ f = ∂f ∂r e r + 1 ∂f ∂θ e θ + ∂f ∂z e z r ∇ · f = 1 ∂rf r + 1 ∂f θ ∂θ + ∂f z r ∂r r ∂z � 1 � � ∂f r � � ∂rf θ � ∂f z ∂θ − ∂f θ ∂z − ∂f z e θ + 1 − ∂f r ∇ × f = e r + e z r ∂z ∂r r ∂r ∂θ ∆ f = ∂ 2 f ∂ 2 f ∂θ 2 + ∂ 2 f ∂r 2 + 1 r 2 ∂z 2
Tensor product � � f k,n x k y n f ( x, y ) = f ( x, y ) = f k,n T k ( x ) T n ( y ) or, better k,n k,n Cylindrical coordinates � f k,m r k e imθ f ( r, θ ) = k,m Despite their seemingly innocuous form, these are not analytic at the origin! f ( r, θ ) = r f ( r, θ ) = cos( θ )
x k y n = ( r cos θ ) k ( r sin θ ) n � � k � � n e iθ + e − iθ e iθ − e − iθ = r k + n 2 2 i � k � r � k + n k � � ( − i ) n e ik ′ θ e − i ( k − k ′ ) θ = k ′ 4 k ′ =0 � n n � � e in ′ θ ( − 1) n − n ′ e − i ( n − n ′ ) θ n ′ n ′ =0 � k � � n � r � k + n k n � � � ( − i ) n ( − 1) n − n ′ e i (2 k ′ − k +2 n ′ − n ) θ = k ′ n ′ 4 k ′ =0 n ′ =0 Wavenumber in θ multiplying r k + n is same parity and restricted to − ( k + n ) ≤ 2 k ′ − k + 2 n ′ − n ≤ k + n j ∞ ∞ ∞ � � � � f jm r j e imθ = f jm r j e imθ f ( r, θ ) = m = −∞ j =0 m = − j j = | m | j + m even j + m even
f = r 2 + 1 f = r 2 − 1 10 r 4 cos(4 θ ) 2 cos(4 θ ) × NO � OK Trigonometric functions with wavenumber m contain m oscillations. As r decreases, oscillations are compressed over decreasing circumference. Requirement that radial function multiplying e imθ begin with r m (that it have an m -th order zero) leads to sufficiently fast damping of oscillation amplitude near origin.
Same result can be demonstrated via differentiation: ∂ ∂x ( x k y n ) = kx k − 1 y n Not singular at x = 0 , since k − 1 < 0 → k = 0 . But: � � ∂ 1 ∂ � r j e imθ � � r j e imθ � ∇ = e r ∂r + e θ ∇ r ∂θ = (e r j + e θ im ) r j − 1 e imθ Singular at r = 0 for j = 0 , m � = 0 . (Require j ≥ m .) � ∂ 2 � � ∂ 2 ∂r 2 + 1 ∂r + 1 ∂ � r j e imθ � r j e imθ � ∆ = r 2 ∂θ 2 r = ( j ( j − 1) + j − m 2 ) r j − 2 e imθ = ( j 2 − m 2 ) r j − 2 e imθ = ⇒ need m = ± j for j = 0 , 1 = ( j 2 − m 2 )(( j − 2) 2 − m 2 ) r j − 4 e imθ ∆ 2 � r j e imθ � Need m = ± j or m = ± ( j − 2) for j ≤ 4 . Etc.
For f to be infinitely differentiable, require j ∞ ∞ ∞ � � � � f jm r j e imθ = f jm r j e imθ f ( r, θ ) = j =0 m = −∞ m = − j j = | m | j + m even j + m even r j � � e imθ � � e ± ijθ + e ± i ( j − 2) θ + . . . r | m | + r | m | +2 + . . .
• We know that monomials r j are a badly conditioned basis. They are small almost everywhere in the interior. Matrix transforming between f ( r j ) and monomial coeffcients is badly conditioned. • What about r m T j ( r ) (called Roberts functions)? ∞ ∞ � � r m f jm T j ( r ) e imθ f ( r, θ ) = m = −∞ j = 0 j + m even Also badly behaved, again because r m exaggerate the boundary. • Bessel functions: � J m � m ∞ � � λr ( ∓ λr/ 2) 2 j � ( λr ) = I m 2 j ! Γ( m + j + 1) j =0 Eigenfunctions of the Laplacian: ∆ J m ( λr ) e imθ = − λ 2 r 2 J m ( λr ) e imθ Correct relations between r exponent m + 2 j and θ wavenumber m . But convergence rate of coefficients in Bessel series is only algebraic.
j ( r ) e imθ = r m P 0 ,m ( j − m ) / 2 (2 r 2 − 1) • One-sided Jacobi basis W m Has good propreties, but too difficult to deal with. • Instead, impose only parity = ⇒ f not analytic: ∞ ∞ ( f, ∇ f regular � � f jm T j ( r ) e imθ f ( r, θ ) = but ∆ f ∼ e 2 iθ /r 2 ) m = −∞ j = 0 j + m even Turns out to be wasteful but not harmful. Coefficient of non-analytic functions are carried around and computed but do not mix with the coefficients of the analytic basis functions. This is even true if parity is not imposed. (Then 3/4 of the functions are non-analytic.)
If using finite differences in r , require at r = 0 ⇒ d f m f m = a + cr 2 + . . . (no linear term) for m even = dr ( r = 0) = 0 for m odd f m = br + . . . = ⇒ f m ( r = 0) = 0 (no constant term) How to incorporate BCs at r = r 0 = 0 into finite difference operator? d 2 f dr 2 ( r 2 ) ≈ αf ( r 3 ) + βf ( r 2 ) + γf ( r 1 ) d 2 f dr 2 ( r 1 ) ≈ αf ( r 2 ) + βf ( r 1 ) f ( r 0 ) = 0 d 2 f f ′ ( r 0 ) = 0 = dr 2 ( r 1 ) ≈ αf ( r 2 ) + ( β + γ ) f ( r 1 ) ⇒ f ( r 0 ) = f ( r 1 ) Can use Cartesian five or nine-point finite-difference stencil at r = 0 , polar coordinates elsewhere. Or omit point at r = 0 .
Full disk Annulus: no singularity BC at r out and regularity at r in Use finite differences or Chebyshev polynomials in r and Fourier series in θ Cluster points at outer boundary, not at the origin Map [ r in , r out ] to [ − 1 , 1] BCs at r in , r out . Either r ∈ [0 , 1] and θ ∈ [0 , 2 π ] or r ∈ [ − 1 , 1] and θ ∈ (0 , π ]
What about vector components ( u r , u θ , u z ) ? All of ( u x , u y , u z ) are like scalars and must obey rules above. u r = cos( θ ) u x + sin( θ ) u y u θ = − sin( θ ) u x + cos( θ ) u y Expanding u x and u y leads to | j +1 | ∞ ∞ ∞ � � � � u r,θ = jm r j e imθ = u r,θ u r,θ jm r j e imθ j =0 m = −∞ m = −| j + 1 | j = | m − 1 | j + m odd j + m odd
Vector Laplacian couples u r and u θ : u r ∆ − 1 − 2 u r g r r 2 ∂ θ 0 → − r 2 ≡ 2 ∆ − 1 = u θ u θ g θ ∆ r 2 ∂ θ 0 r 2 u z u z g z 0 0 ∆ Diagonalize the two-by-two block: � I � u + � u r � � � iI u ± ≡ u r ± iu θ E ≡ = E u − u θ I − iI � ∆ − 1 � r 2 + 2 i E − → r 2 ∂ θ 0 ∆ E − 1 = ∆ − 1 r 2 − 2 i 0 r 2 ∂ θ � u r � g r � � − → ∆ = u θ g θ � u r � g r � � E − → ∆ E − 1 E = E u θ g θ � ∆ − 1 � u + � g + � � � r 2 + 2 i r 2 ∂ θ 0 = ∆ − 1 r 2 − 2 i u − g − 0 r 2 ∂ θ
Spherical Coordinates and Spherical Harmonics Spherical harmonics: Y ℓ,m = Ne imφ P m ℓ (cos θ ) Behavior near poles: ℓ (cos θ ) ∼ sin m θ P m Spherical θ at poles is like r near center of disk: polar cap Many useful relations, such as: ∆ Y ℓ,m = − ℓ ( ℓ + 1) Y ℓ,m r 2 ℓ ( x ) = ( − 1) m (1 − x 2 ) m/ 2 d m P m dx m P 0 ℓ ( x ) (1 − x 2 ) d 1 � � dxP m ( ℓ + 1)( ℓ + m ) P m ℓ − 1 − ℓ ( ℓ − m + 1) P m = ℓ ℓ +1 2 ℓ + 1 � 2 π Y m ℓ ( θ, φ ) Y m ′ δ ℓ,ℓ ′ δ m,m ′ = sin θ dθ dφ ℓ ′ 0
ℓ = 2 ℓ = 4 ℓ = 8 m = 0 m � = 0 , ℓ m = ℓ zonal tesseral sectoral
P m ℓ (cos θ ) ℓ = 5 → m = 1 m = 0 m = 2 As m increases, roots/extrema/variation of polynomials concentrate at ori- gin (equator). Counteracts accumulation of longitude lines ( e imφ ) at poles. ⇒ Areas sampled equally over entire sphere via oscillations of P m = ℓ
Pseudospectral method: transform to grids N ℓ N N � � � � ℓ (cos θ ) e imφ = f m ℓ P m f m ℓ P m e imφ f ( θ, φ ) = ℓ (cos θ ) ℓ =0 m = − ℓ m = − N ℓ = | m | � �� � f m ( θ ) � 2 π N φ � f ( θ, φ ) e − imφ dφ ≈ f ( θ, φ j ) e − imφ j ∆ φ f m ( θ ) = 0 j =1 � π N θ ( m ) � f m f m ( θ ) P m f m ( θ j ) P m 2 πj = ℓ (cos θ ) sin θ dθ ≈ ℓ (cos θ j ) sin θ j ℓ N θ ( m ) 0 � �� � j =1 � FFT in φ direction ∆ θ j Transform to physical grid via weighted sum (matrix mult) in θ direction Grid is equally spaced in φ but more concentrated in θ near equator. ℓ (cos θ ) / sin m θ , Optimal grid for each m would be N θ ( m ) roots of P m but want same θ grid for each m , so use N roots of Legendre poly P 0 N . Retain f m ℓ for ℓ ∈ [ m, N ]
Orientation and role of m ℓ = 0 ℓ = 1 → ℓ = 2 ℓ = 3 ℓ = 4 = − Ne ± iφ sin θ = N ∓ x − iy = N x � � Y ± 1 Y − 1 1 − Y 1 1 1 1 2 r r √ √ 2 z = N ± z − iy � � Y − 1 Y 0 ± 1 2 Y 0 1 + 1 + Y 1 1 = N 2 cos θ = N √ 2 1 1 r r So m depends on orientation w.r.t. z -axis as well as on variation, unlike ℓ . = − ℓ ( ℓ +1) Rotating ( z, x ) → ( − x, z ) changes m , but ∆ Y m Y m ∀ m ℓ ℓ r Distinguished choice of z axis if there is rotation.
Fornberg: finite differences on an equally spaced grid Pole � � ∂ 2 f 1 ∂ θ∂f 1 cos ˜ ∆ surface f = + cos ˜ ∂ ˜ ∂ ˜ cos ˜ ∂φ 2 θ θ θ θ Here ˜ θ is latitude, measured from equator h ≡ ∆ φ , k ≡ ∆˜ θ
Recommend
More recommend
