Spectral solver for Cauchy problems in polar coordinates using discrete Hankel transforms Rundong Zhou13and Nicolas Grisouard2

2025-05-03 0 0 2.27MB 24 页 10玖币

侵权投诉

Spectral solver for Cauchy problems in polar

coordinates using discrete Hankel transforms

Rundong Zhou1,3* and Nicolas Grisouard2

1*Division of Engineering Science, University of Toronto, 40 St. George

Street, Toronto, M5S 2E4, Ontario, Canada.

2Department of Physics, University of Toronto, 60 St. George Street,

Toronto, M5S 1A7, Ontario, Canada.

3Now at Department of Physics, Chalmers University of Technology,

Kemigården 1, Göteborg, SE-412 96, Sweden.

*Corresponding author(s). E-mail(s): rundongz@student.chalmers.se;

Contributing authors: nicolas.grisouard@utoronto.ca;

Abstract

We introduce a Fourier-Bessel-based spectral solver for Cauchy problems fea-

turing Laplacians in polar coordinates under homogeneous Dirichlet boundary

conditions. We use FFTs in the azimuthal direction to isolate angular modes,

then perform discrete Hankel transform (DHT) on each mode along the radial

direction to obtain spectral coeﬃcients. The two transforms are connected via

numerical and cardinal interpolations. We analyze the boundary-dependent error

bound of DHT; the worst case is ∼N−3/2, which governs the method, and

the best ∼e−N, which then the numerical interpolation governs. The complex-

ity is O[N3]. Taking advantage of Bessel functions being the eigenfunctions of

the Laplacian operator, we solve linear equations for all times. For non-linear

equations, we use a time-splitting method to integrate the solutions. We show

examples and validate the method on the two-dimensional wave equation, which

is linear, and on two non-linear problems: a time-dependent Poiseuille ﬂow and

the ﬂow of a Bose-Einstein condensate on a disk.

Keywords: Spectral methods, Discrete Hankel transforms, Initial value problems,

Nonlinear partial diﬀerential equations, Bose-Einstein condensates

arXiv:2210.09736v4 [physics.comp-ph] 23 Jul 2023

1 Introduction

The Laplacian operator is associated with many important physical problems, includ-

ing diﬀusion, the Schrödinger equation or the wave equation. It is more complicated

under polar and spherical coordinates than in Cartesian coordinates. The literature

has explored the use of the ﬁnite diﬀerence method for solving Poisson-type equations

in cylindrical and spherical geometries [1,2], and spectral approaches are often used

to avoid coordinate singularities [3,4]. In this paper, we introduce a novel spectral

solver suited for time-dependent problems featuring Laplacian in polar coordinates.

We consider a function ψon a unit disk r∈[0,1] that satisﬁes the homogeneous

Dirichlet boundary condition ψ(1, θ)=0. The periodicity in the azimuthal direction

allows for a decomposition of ψas a Fourier series, namely,

ψ(r, θ) =

∞

q=−∞

fq(r)eiqθ,with fq(r) = 1

2πZ2π

ψ(r, θ)e−iqθ dθ(1)

and where qdenotes the angular mode number. There are various choices of basis func-

tions to decompose the radial function fq(r). Boyd and Yu [5] list ﬁve spectral methods

that use decompositions of the radial function. Among the listed basis functions, four

are families of orthogonal polynomials, such as Zernike [6], Logan-Shepp [7], modi-

ﬁed Chebyshev [3,8,9] or modiﬁed Jacobi polynomials [10]. Most of them focus on

solving boundary value problems and time-independent partial diﬀerential equations

(PDEs) in polar coordinates, while Cauchy problems involving time evolution receive

comparatively little attention.

To build a Laplacian solver suitable for time-dependent problems, instead of choos-

ing orthogonal polynomial bases, we decided to use Bessel functions of the ﬁrst kind as

our basis to decompose the radial function fq(r), given by the Fourier-Bessel series [11],

fq(r) =

∞

j=1

aq,j Jq(kq,j r),with aq,j =2

q+1(kq,j )Z1

rfq(r)Jq(kq,j r)dr. (2)

Jqis the qth-order Bessel function of the ﬁrst kind and kq,j denotes its jth non-negative

zero. Temme [12] provides a fast and accurate algorithm to compute these roots.

The Fourier-Bessel basis (i.e., the combination of the azimuthal Fourier basis of

Eqn. 1and radial Bessel basis of Eqn. 2) has an algebraic rate of convergence ≲1/N5/2

for functions satisfying homogeneous Dirichlet boundary conditions [5], where Nis

the number of basis functions used for approximation. This compares unfavourably

with many orthogonal polynomial bases which have an exponential rate. However,

Fourier-Bessel modes are eigenfunctions of the Laplacian operator, and the boundary

conditions are enforced by the basis functions themselves. Such virtues ensure that

a pure spectral scheme can be applied at each iteration. It makes the Fourier-Bessel

basis a competitive choice for solving time-dependent initial value problems associated

with Laplacians under homogeneous Dirichlet conditions in polar coordinates. We

further notice a recent related work [13] on obtaining Fourier-Bessel coeﬃcients but

considering time-independent problems on a Cartesian sampling grid.

To decompose the radial function fq(r)into Fourier-Bessel series, we use the

discrete Hankel transform (DHT). It was ﬁrst introduced mathematically by John-

son [14], then independently re-invented by Yu et al. [15] for the zeroth order mode

(q= 0) and Guizar-Sicairos et al. [16] for all integer orders. Finally, It was categorized

by Baddour [17] as a discrete variation of general Fourier transform.

This paper is organized as follows: we ﬁrst reformulate the DHTs as pseudospectral

collocation method by introducing their discrete inner products, quadrature weights,

pseudospectral grid points, and cardinal interpolations in Section 2.1 and 2.2. We

analyze error and complexity and discuss the boundary dependency of the error bound

in Section 2.3. Then we introduce a systematic way to apply the method to compute

Laplacian and solve nonlinear time-dependent equations in Section 2.4. In Section 3

we show three examples: the linear 2-D wave equation where we test the convergence

of the method numerically (§ 3.1), a time-dependent Poiseuille ﬂow equation (§ 3.2),

and the Gross-Pitaevskii equation (§ 3.3). In Section 4we conclude the paper and

discuss potential improvements.

2 Method

2.1 Discrete Hankel transform and pseudospectral grid points

Our ﬁrst step is to reformulate the discrete Hankel transforms as pseudospectral collo-

cation methods. The optimal pseudospectral grid points are the zeros of the (N+ 1)th

basis function ϕN+1(r)[18, § 4.3], where Nis the total number of basis functions, or

radial modes, used to approximate the radial function fq(r). The basis functions of a

qth-order Hankel transform are then

ϕq,i(r) = Jq(kq,ir).(3)

Thus, the (N+ 1)th basis function is Jq(kq,N+1r), and we choose its pseudospectral

grid points to be

{rq,i}i=1,2,...,N =kq,i

kq,N+1

.(4)

We can now deﬁne the discrete inner product as [19, § 3.1.4]

(f, g)q≡

i=1

wq,if(rq,i)g(rq,i),where wq,i =2

q,N+1J2

q+1(kq,i)(5)

are the quadrature weights.

Bessel functions are orthogonal under our discrete inner product, namely [14]

q+1(kq,m)Jq(kq,mr), Jq(kq,nr)q=δmn.(6)

The discrete Hankel transform of the radial function fq(r)is now deﬁned through the

discrete inner product as

Fq,j =fq(r), Jq(kq,j r)q.(7)

By substituting the deﬁnition of discrete inner product and pseudospectral grid points

{rq,i}, the forward discrete Hankel transform formula is then given by

Fq,j =2

q,N+1

i=1

Jq(kq,ikq,j /kq,N+1)

q+1 (kq,i)fq(rq,i).(8)

The value of the radial function fq(r)at each pseudospectral grid point rq,i can be

recovered by the Fourier-Bessel series Eqn. (2), namely,

fq(rq,i)=2

j=1

Jq(kq,ikq,j /kq,N+1)

q+1 (kq,j )Fq,j .(9)

The transform formulae Eqns. (8) and (9) demonstrate symmetry and can be

implemented as matrix multiplications. Introducing

fq,i ≡fq(rq,i)and Mq,ij ≡ϕq,i(rq,j )wq,i =2Jq(kq,ikq,j /kq,N +1)

q,N+1J2

q+1 (kq,i),(10)

we obtain ⃗

Fq=Mq·⃗

fq;⃗

fq=k2

N+1Mq·⃗

Fq,(11)

where Mqis unitary, i.e., MqM−1

q=I. The discrete Hankel transform can be catego-

rized as a Matrix Multiplication Transformation [18, § 10.4], and it demonstrates the

spectral methods’ virtue of one-to-one correspondence between spectral coeﬃcients

Fq,j and pseudospectral grid points {fq(rq,i)}.

2.2 Interpolation and cardinal functions

To obtain the spectral coeﬃcients aq,j of Eqn. (2), we perform a Fast Fourier Trans-

form (FFT) along the azimuthal direction followed in the radial direction by one DHT

for each angular mode q. The pseudospectral grid points {rq,i}in Eqn. (4) deﬁned for

the DHTs are not evenly spaced and vary with the angular mode q. Instead of using

Baddour and Yao’s unevenly sampled grid [20,21], which causes unwanted artifacts

in the center of the domain, we choose an evenly sampled grid that is ﬁner in the

radial direction before performing the azimuthal FFT (hereafter referred to as the

FFT grid), as visualized in Figure 1a. We then numerically interpolate each radial

function fqon the corresponding DHT grid {rq,i}i=1,...,N to obtain {fq(rq,i)}i=1,...,N ,

as shown in Figure 1b.

To achieve optimal numerical accuracy for a given N, the number of radial sam-

pling points of the ﬁner FFT grid should be at least 2N(see § 2.3.1) The number of

(a) FFT grid (b) DHT grid

Fig. 1: Diﬀerence between the evenly sampled FFT grid on ψ(r, θ)and the

pseudospectral {rq,i}grid for the DHTs in the case of N= 8. In panel (b), the

sampling points for q= 0 (marked as purple squares) do not contain the origin since

r0,1>0.

azimuthal sampling points is twice that of the desired angular modes q. Conversely,

when we transform from the spectral domain aq,j back to the physical domain ψ(r, θ),

we perform inverse DHTs followed by an inverse FFT. We use the cardinal functions

of Fourier-Bessel series to interpolate from the pseudospectral grid {rq,i}to the ﬁner

FFT grid [17], namely,

fq(r) =

∞

i=1

fq(rq,i)Cq,i(r),where Cq,i(r) = 2kq,i

q,i −r2k2

q,N+1

Jq(rkq,N+1)

Jq+1(kq,i).(12)

We discuss the error introduced by these two interpolations in the following section.

2.3 Error analysis and complexity

2.3.1 Numerical interpolation

The numerical interpolation connecting the ﬁner FFT grid to the pseudospectral grid

{rq,i}introduces an error that largely depends on the resolution of the ﬁner FFT grid

and the choice of numerical method. We choose the cubic spline interpolation, whose

global error EIto approximate the radial function fq(r)is bounded by [22]

EI≤5

384 f(4)

qh4,

where his the size of the radial interval of the FFT grid. According to our sam-

pling strategy explained in § 2.2, there are 2Nradial FFT sampling points and hence

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

SpectralsolverforCauchyproblemsinpolarcoordinatesusingdiscreteHankeltransformsRundongZhou1,3*andNicolasGrisouard21*DivisionofEngineeringScience,UniversityofToronto,40St.GeorgeStreet,Toronto,M5S2E4,Ontario,Canada.2DepartmentofPhysics,UniversityofToronto,60St.GeorgeStreet,Toronto,M5S1A7,Ontario,Canada...

展开>> 收起<<

Spectral solver for Cauchy problems in polar coordinates using discrete Hankel transforms Rundong Zhou13and Nicolas Grisouard2.pdf

共24页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

Spectral solver for Cauchy problems in polar coordinates using discrete Hankel transforms Rundong Zhou13and Nicolas Grisouard2

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: