Numerical Solution of an Extra-wide Angle Parabolic Equation through Diagonalization of a 1-D Indefinite Schrödinger Operator with a Piecewise Constant Potential

2025-05-02 0 0 1.41MB 26 页 10玖币

侵权投诉

Numerical Solution of an Extra-wide Angle Parabolic

Equation through Diagonalization of a 1-D Indeﬁnite

Schr¨odinger Operator with a Piecewise Constant

Potential

Sarah D. Wrighta,∗, James V. Lambersa

aSchool of Mathematics and Natural Sciences, The University of Southern Mississippi, 118

College Dr #5043, Hattiesburg, MS 39406, USA

Abstract

We present a numerical method for computing the solution of a partial diﬀer-

ential equation (PDE) for modeling acoustic pressure, known as an extra-wide

angle parabolic equation, that features the square root of a diﬀerential operator.

The diﬀerential operator is the negative of an indeﬁnite Schr¨odinger operator

with a piecewise constant potential. This work primarily deals with the 3-piece

case; however, a generalization is made the case of an arbitrary number of

pieces. Through restriction to a judiciously chosen lower-dimensional subspace,

approximate eigenfunctions are used to obtain estimates for the eigenvalues of

the operator. Then, the estimated eigenvalues are used as initial guesses for the

Secant Method to ﬁnd the exact eigenvalues, up to roundoﬀ error. An eigenfunc-

tion expansion of the solution is then constructed. The computational expense

of obtaining each eigenpair is independent of the grid size. The accuracy, eﬃ-

ciency, and scalability of this method is shown through numerical experiments

and comparisons with other methods.

Keywords:

Partial diﬀerential equations, Schr¨odinger operator, Secant Method,

Extra-wide angle Parabolic equation, Wave propagation

1. Introduction

The problem to be considered in this paper is an extra-wide angle parabolic

equation (termed EWAPE2 in [1]), a PDE that is applicable to sound wave

propagation where backscattered energy is negligible [2]. The PDE is expressed

in cylindrical coordinates as

dr =±i r∂2

∂z2+α2I!u0< z < π, r > 0,(1)

∗Corresponding author.

Preprint submitted to Appl. Numerical Math. October 10, 2022

arXiv:2210.03655v1 [math.NA] 7 Oct 2022

where α(r, z) is deﬁned to be

α=ω

c(r, z).(2)

Here, ωis the angular frequency, and c(r, z) is the sound speed as a function of

the range, r, and depth, z. The range will be treated as a temporal variable for

the purpose of numerical solution.

This PDE has many applications in physics [1]; our interest is in modeling

acoustic pressure in the ocean, in the case where the sound speed is weakly

dependent on the range r. In this context, the sound speed has relatively small

variation from its average with respect to depth [3]. Therefore, we treat the

sound speed as constant within each of two boundary layers, and within a third

interior region. That is, we treat the wavenumber as a piecewise constant coef-

ﬁcient, α(z), over three pieces. Previous work has demonstrated that solution

methods using eigenfunction expansion have been eﬀective in solving PDEs with

piecewise constant coeﬃcients [4]. In this paper, highly accurate eigenvalues and

eigenfunctions will be computed, and then used to construct the solution of our

version of the EWAPE2, in a way that sidesteps the diﬃculties encountered

by numerical methods due to the need to approximate the square root of a

diﬀerential operator. Although this work is connected to previous work on re-

lated eigenvalue problems [4, 5], a key diﬀerence is that two distinct approaches

are required to obtain reasonably accurate initial guesses for the eigenvalues,

that can then be improved via iteration–one for low-frequency eigenfunctions,

and one for high-frequency. By exploiting the orthogonality of eigenspaces, an

eigenfunction of either kind can be obtained in O(1) ﬂoating-point operations.

The outline of the paper is as follows. In Section 2, we will discuss the

background of the problem in further detail. Section 3 will describe the approach

to ﬁnding eigenfunctions. Section 4 will discuss the solution of the PDE through

an eigenfunction expansion. Section 5 is a compilation of all of the numerical

results obtained using our algorithm. This will include evaluation of accuracy

and eﬃciency and comparison of our results to those of other numerical methods.

Section 6 will describe a generalization of our algorithm to an arbitrary piecewise

constant wavenumber. Section 7 will conclude the paper with ﬁnal thoughts

regarding the proposed algorithm, and possible future work.

2. Background

The Helmholtz equation is an eigenvalue problem for the Laplacian oper-

ator and has many applications in physics such as seismology, acoustics, and

electromagnetic radiation [6]. In cylindrical coordinates, it has the form

∂2

∂r2+1

∂

∂r +∂2

∂z2+α2P= 0.

This equation can be factored into two equations that model either forward

propagation or backward propagation. The resulting PDE, also known as the

“one-way equation”, is deﬁned as

dr =±i√Lu, 0< z < π, r > 0,(3)

where the spatial domain (that is, the range of depths) is chosen to be (0, π)

solely for notational convenience, and the diﬀerential operator Lis deﬁned by

L=∂2

∂z2+α2I. (4)

This one-way equation only models waves traveling in one direction and so we

must assume that any propagation in the opposite direction, such as echoes, is

negligible.

We also impose the “initial” condition

u(z, 0) = f(z),0< z < π, (5)

and homogeneous Dirichlet boundary conditions. The solution uto this one-

way equation is related to the solution, P, of the Helmholtz equation by the

change of variable

u=√rP.

We are assuming that the sound speed depends only on the depth and not

range, to justify the factoring of the Helmholtz equation [2]. Because the sound

speed ﬂucutates very little relative to its average, we can approach this problem

as if it were nearly piecewise constant. Therefore, for our simulation, we will

model α2as a piecewise constant function consisting of three pieces, where the

ﬁrst and last piece come from the wavenumber at the bounding depths. The

average of the wavenumber over the interior depths will be used for the middle

piece.

The one-way equation (3) includes the square root of a negated 1-D Schr¨odinger

operator. A general Schr¨odinger operator is the sum of the kinetic energies, rep-

resented as a negated Laplacian, and the potential energies, accounted for by V

[7]:

H=−



j=1

∂2

∂x2

j

+V.

Our operator Lin the one-way equation is a scalar multiple of a Schr¨odinger

operator with a negative potential; such Schr¨odinger operators were studied

in [8]. Taking the square root of Lposes substantial diﬃculties for numerical

methods. For a narrow-angle parabolic equation, √Lis approximated by using

a ﬁrst-order Taylor expansion of

α0√1 + `

where

`=L

α2

0−1.

This was originally proposed by Claerbout for extrapolation of seismic waves

[9]. In order to obtain a suﬃciently accurate approximation, Taylor expansion

is used, which yields √1 + `≈1 + `

By limiting the Taylor expansion to only the ﬁrst two terms, this approach sac-

riﬁces accuracy for eﬃciency. This approach will not work well for our problem

because the matrix Aobtained by discretizing Lis ill-conditioned. For a wide-

angle parabolic equation, a rational approximation is used instead [2], but the

applicability of this approach is also limited, again due to the ill-conditioning

of this matrix.

Another approach is to use the Schur decomposition A=QDQH, where Qis

unitary, and Dis diagonal, as Ais symmetric. It follows that √A=Q√DQH.

This method is very expensive, as it requires approximately 281

3n3ﬂops [10],

where Ais n×n. A far more eﬃcient approach would be to compute accurate

eigenvalues and eigenfunctions of the operator directly, without relying on any

spatial discretization. We will develop such an approach in this paper.

2.1. Previous Work

Similar to previous work concerning eigenfunction expansion [4, 5], the deriva-

tion of the proposed algorithm will begin by deﬁning the eigenfunctions as com-

binations of sine and cosine functions. Each eigenfunction must satisfy bound-

ary conditions, and continuity requirements at the interfaces between pieces.

The system of equations arising from these conditions can then be reduced by

eliminating most parameters. Once the system is reduced to a single nonlinear

equation, we can solve for the remaining parameter using the Secant Method

[11]. This remaining variable is the frequency corresponding to the third piece

of the eigenfunction. In order to use The Secant Method, we will seek sim-

ple approximations of each eigenvalue, to serve as initial guesses, as was done

in [4, 5], but as the diﬀerential operator is of a diﬀerent form, a substantially

diﬀerent approach to obtaining these initial guesses will be required.

2.2. Behavior of Eigenfunctions

To guide the development of our algorithm, we use Matlab to compute

eigenvalues and eigenvectors of a matrix that discretizes our operator L, so that

we may understand their behavior. Figures 1, 2, and 3 show the graphs of

the absolute value of the smallest several eigenvalues. The ﬁrst few are positive

eigenvalues that are bounded in magnitude; however, the remaining eigenvalues

are negative and are unbounded.

Each plot within Figure 4 displays an eigenfunction in blue and the function

sin(jz), for some j, in red. We can see that eigenfunctions that correspond

to larger eigenvalues can be approximated by a sine function; however, as the

eigenvalues approach 0, they become more diﬃcult to approximate using such

a sine function. In fact, when the eigenvalues become positive, they nowhere

resemble a sine function, as we can see from Figure 5, so we must therefore

develop a diﬀerent approach for approximating these eigenfunctions.

Figure 1: Absolute values of all eigenvalues of “one way” equation with α1= [2,1,2], ρ=

[1/3,2/3], and N= 2048. The positive eigenvalues are seen in red and the eigenvalues that

correspond to negative eigenvalues are blue.

3. Methodology

The eigenfunctions for the piecewise constant coeﬃcient case, with three

pieces, will have the form

Vj(z) = 





Vj1(z) = Aj1cos(ωj1z) + Bj1sin(ωj1z) 0 ≤z < πρ1,

Vj2(z) = Aj2cos(ωj2z) + Bj2sin(ωj2z)πρ1≤z < πρ2,

Vj3(z) = Aj3cos(ωj3z) + Bj3sin(ωj3z)πρ2≤z < π,

where ωj1, ωj2, ωj3are the unknown frequencies and Aj1, Aj2, Aj3and Bj1, Bj2, Bj3

are related to the unknown amplitudes and phase shifts. A system of linear

equations that characterizes the eigenfunction comes from imposing conditions

on Vj(z). These conditions are homogeneous Dirichlet boundary conditions and

continuity conditions. The continuity conditions require that the one-sided lim-

its of eigenfunctions at the interfaces are equal, as well as those of their ﬁrst

derivatives [12]. We will also set conditions for how the frequencies ωj1,ωj2,

and ωj3relate to each other.

•Dirichlet Boundary Conditions:

Vj1(0) = 0,(6)

Vjn(π)=0.(7)

•Continuity:

Vj1(πρ1) = Vj2(πρ1),(8)

j1(πρ1) = V0

j2(πρ1),(9)

Vj2(πρ2) = Vj3(πρ2),(10)

j2(πρ2) = V0

j3(πρ2).(11)

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

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

10 玖币 0人已下载

立即下载

摘要：

NumericalSolutionofanExtra-wideAngleParabolicEquationthroughDiagonalizationofa1-DIndeniteSchrodingerOperatorwithaPiecewiseConstantPotentialSarahD.Wrighta,,JamesV.LambersaaSchoolofMathematicsandNaturalSciences,TheUniversityofSouthernMississippi,118CollegeDr#5043,Hattiesburg,MS39406,USAAbstractWepr...

展开>> 收起<<

Numerical Solution of an Extra-wide Angle Parabolic Equation through Diagonalization of a 1-D Indefinite Schrödinger Operator with a Piecewise Constant Potential.pdf

共26页,预览5页

还剩页未读，继续阅读

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

Numerical Solution of an Extra-wide Angle Parabolic Equation through Diagonalization of a 1-D Indefinite Schrödinger Operator with a Piecewise Constant Potential

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: