高阶微分方程

2024-05-17

高阶微分方程（共6篇）

高阶微分方程篇1

本文研究如下的高阶微分方程

undefined

的周期解的存在性,(1)式中g,p都是定义在R上的实连续函数,b≠0,τ≥0,p以T为周期,且

∫undefinedp(x)dx=0。本文利用重合度理论获得了(1)式至少存在一个T-周期解的充分条件,其结果是如下定理.

定理如果下列条件成立:

(ⅰ)存在正常数M,使得g(x)≤M,∀x∈R;

(ⅱ)undefined;

则方程(1)至少存在一个T(T>0)周期解。

证明考察方程

undefined

这里λ∈(0,1),设x(t)是方程(2)的任一T-周期解,将方程(2)两边同时从0到T积分得

undefined

因此存在t0∈(0,T),使得 bx(t0)+g(x(t0-τ))=0,从而有

undefined

于是

undefined

因为

undefined

k为正整数,

于是方程(2)两边同乘以x(t),再从0到T积分得

undefined

因此

undefined

(4)式中undefined|。

因为x(0)=x(T),则存在t1∈[0,T],使得x′(t1)=0,有

undefined

同理x′(0)=x′(T),存在t2∈[0,T],使得x″(t2)=0,有

undefined

从而∫undefined|x′(t)|dt≤T∫undefined|x″(t)|dt≤…≤Tn-1×

∫undefined|x(n)(t)|dt。所以不等式(4)成为

undefined

因此由条件(ⅱ)及上式知,存在与λ无关的数R1>0,使得

undefined

进而可得存在与λ无关的数R2>0,使得

undefined

从而由方程(2)有

undefined

由(5)式及(7)式知必存在与λ无关的数rj>0(j=1,2,…,2n-1),使得

undefined

取undefined,令X={x(t)∈C(R,R)|x(t+T)=x(t)},undefined,Nx(t)=-cx′-bx-g(x(t-τ))+p(t)。这时有Ker L=R,同时定义投影算子为

P:undefined。Q:X→X/lmundefined。

则Ker L=lm P,Ker Q=lm L,即L是指标为零的Fredholm算子,且可证明N在undefined⊂X上L-紧。方程(2)即为算子方程Lx=λNx,λ∈(0,1)。根据对(2)式周期解界的估计及已知条件,有

1)Lx≠λNx,∀x∈Ker L∩∂Ω,还需证明

2) QNx≠0,∀x∈Ker L∩∂Ω;

3) deg(QN,Ω∩Ker L,0)≠0。

事实上,当x∈Ker L∩∂Ω时,x为常数且

undefined,有

undefined

即QNx≠0,于是2)成立;作变换:F(x,μ)=μbx+(1-μ)[bx+g(x(t-τ))],对任意x∈Ker L∩∂Ω,μ∈[0,1],x为常数且undefined。

undefined

。所以当undefined时,F(x,μ)≠0,同理可证当undefined时,F(x,μ)≠0,因而F(x,μ)为同伦变换,因此

deg{QN,Ω∩Ker L,0}=deg{-bx-g(x(t-τ)),Ω∩Ker L,0}=deg{-bx,Ω∩Ker L,0}≠0。

故3)成立,由重合度理论知,方程(1)至少有一个T(T>0)周期解。

文献[1]的结果是本文定理的简单推论。事实上,只要在本文的方程(1)中取n=1即明。在本文的方程(1)中取n=1,c=0,即可得到文献[2]的结果。

摘要：考虑一类高阶微分方程ax(2n)(t)+cx′(t)+bx(t)+g[x(t-τ)]=p(t),利用重合度理论,获得了此类方程至少存在一个T-周期解的充分条件。

关键词：高阶微分方程,周期解,重合度

参考文献

[1]唐美兰,刘心歌,刘心笔,等.一类时滞Duffing型方程周期解的存在性.上饶师范学院学报,2004;24(6):12—15

[2]张正球,庾建设.一类时滞Duffing型方程周期解.高校应用数学学报,1998;13(4):389—392

高阶微分方程篇2

一类亚纯函数系数的高阶非齐次线性微分方程解的增长性

主要研究了一类亚纯函数系数的.高阶非齐次线性微分方程无穷级亚纯解的增长性问题,对大多数亚纯解的超级、二级不同零点收敛指数得到了精确估计.

作者：郑秀敏陈宗煊曹廷彬涂金作者单位：江西师范大学,数学与信息科学学院,江西,南昌,330027刊名：江西师范大学学报(自然科学版) ISTIC PKU英文刊名：JOURNAL OF JIANGXI NORMAL UNIVERSITY(NATURAL SCIENCES EDITION)年，卷(期)：28(5)分类号：O174.52关键词：线性微分方程亚纯函数超级二级收敛指数

高阶微分方程篇3

若参数方程

确定y与x的函数关系,此处假设φ( t) ,ψ( t) 都可导且φ'( t) ≠0.

1. 由导数的几何意义推导

P0(φ( t0) ,ψ( t0)) 点在曲线上,切线P0T的斜率可由割线P0P的斜率取极限而得,而割线为

于是切线斜率为

2. 由复合函数与反函数的导数推导

在( 1) 式中,如果函数x =φ( t) 具有单调连续反函数t =φ- 1( x),将其代入y = φ( t) 中,得复合函数y = ψ (φ- 1( x)),则由复合函数与反函数的求导法则,就有

由1,2知,由参数方程所确定的函数y对自变量x的导数,等于函数y对参数t的导数除以自变量x对参数t的导数.

3. 由导数也称微商推导而得

dy = ψ' ( t) dt,dx = φ' ( t) dt,从而得 ,即dy/dx=ψ'( t)/φ'( t).

由参数方程所确定的函数y对自变量x的导数,等于函数y对参数t的微分除以自变量x对参数t的微分.

4. 将 y 视为自变量 t 的函数,x 为中间变量进行推导

二、高阶导数

在一阶导数的推导中,上述四种方法各有优点和不足,第一种方法,几何意义明显,但推导较麻烦,后三种推导简单,但看不出几何意义,一般教材采用前两种方法. 而高阶导数是在一阶导数的基础上定义的,所以上述后三种方法给求高阶导数提供了方法. 我们知道,函数y对自变量x的导数仍是x的函数,所以,由参数方程求得dy/dx=ψ'( t)/φ'( t),从而得到

一阶导数的参数方程为

由参数方程所确定的函数对自变量的导数,等于函数对参数t的导数除以自变量对参数t的导数或由参数方程所确定的函数对自变量的导数,等于函数对参数t的微分除以自变量对参数t的微分.

我们用同样的方法可以继续求d3y/dx3.

一般地,若用上述方法已经求得 ,则 n - 1阶导数的参数方程为

由参数方程所确定的函数对自变量的导数,等于函数对参数t的导数除以自变量对参数t的导数或由参数方程所确定的函数对自变量的导数,等于函数对参数t的微分除以自变量对参数t的微分从而有dny/dxn=Ψ'( t)/φ'( t).

例求下列参数方程所确定的函数的三阶导数

一阶导数的第四种推导,可以作为是参数方程所确定的函数的高阶导数的公式法.

一般地,若dn - 1y/dxn - 1已经求出,则

高阶微分方程篇4

6.1 多元函数微分的基本概念

6.1.3 偏导数6.1.4 高阶偏导数(导学)

一、一元函数导数相关知识

1．某商店卖两种牌子的果汁，本地牌子每瓶进价1元，外地牌子每瓶进价1.2元，店主估计，如果本地牌子的每瓶卖x元，外地牌子的每瓶卖y元，则每天可卖出705x4y瓶本地牌子的果汁，806x7y瓶外地牌子的果汁,问：(1)店主每天的收益为多少？（2）收益对不同价格x，y的变化率为多少？

二、多元函数有关问题

1．偏导数符号“”怎么读？

2．多元函数的偏导数几何意义

3．怎样求偏导数?

4．fx(x,y)与正fx(x0,y0)两者是怎样的关系?

三、举例与练习

1．求ux2y2xy的偏导数。z

2．求函数zx23xy2y2在点(2,1)处的两个偏导数

3．设zxy(x0)，求证xz1z2z yxlnxy

u2zu)()2()21 xyz4．设ux2y2z2，求证（xy22,xy0，求f(0,0)和f(0,0)5．函数f(x,y)xyxy22xy00,6．求函数zx3y3x2y3的二阶偏导数.四、思考题

1.二元函数f(x,y))在点(x0,y0)处的偏导数fx(x0,y0)与一元函数(x)f(x,y0)在点x0处的导数(x0)是否相同?

高阶微分方程篇5

1 一阶声波波动方程

假设体力为零, 三维一阶声波方程有具体形式如 (1) 式所示:

其中, Vx、Vy、Vz为速度分量, P为地压分量, 为介质密度, vÁ为体积模量, v为纵波速度。

2 高阶有限差分格式

在用交错网格有限差分数值模拟时, 将应力和速度分量分别设置在和t时刻进行计算[5]采用的空间网格节点分布如图1及表1所示。

相对于空间高阶差分, 时间上的高阶差分会极大地增加计算量, 而时间频散在地震勘探所需的时间采样上并不明显[6], 因此基于运算效率和精度的综合考虑, 本文在时间上采用二阶差分格式, 空间上采用高阶差分格式。

对于微分算子其2N阶有限差分格式如下[4]:

式中为差分方程待定系数, 它可通过方程 (3) 进行求解。

综合以上分析可得三维一阶声波方程的2N阶交错网格差分格式如下:

3 吸收边界条件

有限差分数值模拟中常采用的完美匹配层 (PML) 边界吸收条件, 其基本思想是将波场分裂成垂直和平行于传播方向的两组, 在人工边界处快速衰平行界面法向传播的平面波, 对其它方向不衰减, 从而达到减小边界反射的目的, 该方法可很好地吸收不同方向不同频率的边界反射[7,8,9]。对于一阶声波方程, 由于速度分量的计算仅是声压分量与其对应的坐标轴的一阶空间导数, 因此无须对其进行分裂处理。只要对声压分量进行分裂, 将其分解为Px, Py, Pz三个部分, 即:

与其对应的PML边界条件方程可写为:

此外, 三个速度位移方程的PML形式:

(6) ~ (11) 式中d (x) 、d (y) 、d (z) 分别为x, y, z三个方向的衰减因子, Collino推导了其具体表达式如下[10]:

式中vmax为最大纵波速度, 为匹配层厚度, R为理想状态下介质的反射系数。综上分析可得三维一阶声波方程的PML边界条件的差分格式如下:

4 稳定性条件

通过傅里叶分析方法, 可得到三维交错网格高阶有限差分格式的稳定性条件如 (5) 式所示[11]。

式中, t为时间采样间隔, Vmax为计算区域内最大纵波速度, x、y、z分别为x、y、z三个方向的网格间距。

5 模型试算

为了验证算法的正确性, 本文采用均匀介质模型进行了数值模拟和分析。图2所示为一均匀介质模型, 模型尺度为1200m×1200m×1200m, 网格剖分大小为10m×10m×10m, 介质速度为2000m/s, 震源位于坐标 (600m, 600m, 600m) 处。采用主频为30Hz的雷克子波用为震源, 时间步长为1ms, 模拟记录长度为1s, 差分阶数为时间2阶, 空间10阶。

图3是未加吸收边界进t=400ms的波场快照, 可见未加边界处理时, 地震波在模型边界处理发生了反射, 干扰了地震记为, 图4是经PML边界处理后地波场快照, 可见经边界处理后, 在边界处, 地震波得到了有效地吸收, 未发生边界反射。

图4是PML边界处理前后半空间均匀介质的地震记录对比。通过对比可以看出边界处理前, 直达波在边界处理发生了多次反射, 干扰了地震记录, 而经过边界处理后, 地震记录中的边界反射得到了有效吸收, 且不受任何角度和频率的限制。

6 结论

本文从理论上推导了三维一阶声波方程高阶交错网格有限差分正演的差分格式, 并采用完全匹配边界条件对模型边界进行了处理。正演结果表明, 该方法具有空间频散小的特点, 可以很好地吸收衰减模拟边界产生的边界反射, 不受反射角度和频率范围的限制。

摘要：在地震波数值模拟过程中, 采用高阶交错网格有限差分法可以有效地压制频散, 提高模拟精度。本文推导了三维一阶声波方程的空间任意偶数阶有限差分格式, 并利用完全匹配边界条件对三维地质模型边界进行处理。数值模拟结果表明, 该方法可以有效地吸收人为边界反射, 不产生任何边界反射, 取得了良好的吸收效果。

关键词：一阶声波方程,交错网格,PML边界条件

参考文献

[1]杜增利, 李亚林, 尹成, 高宏亮.虚谱法一阶应力-速度方程地震模拟[J].石油地球物理勘探, 2009.10, 44 (5) :637.

[2]陈耿毅, 余钦范, 蔡希玲等.地震模拟技术新进展-第67届EAGE年会论文综述[J].勘探地球物理进展, 2005, 28 (6) :439~448.

[3]张永刚.地震波场数值模拟方法[J].石油物探, 2003.6, 42 (2) :143~148.

[4]和少伟.弹性介质地震波场的数值模拟[D].荆州:长江大学.2012.4:1~2.

[5]董良国, 马在田, 曹景忠, 王忠华等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报, 2000.5, 43 (3) :413.

[6]邓帅奇.全空间弹性波场数值模拟与逆时偏移成像方法研究[D].徐州:中国矿业大学.2012.5:23.

[7]牟永光, 裴正林.三维复杂介质地震数值模拟[M].北京:石油工业出版社, 2005:182.

[8]王永刚, 刑文军, 谢万学, 朱兆林.完全匹配层收边界条件的研究[J].中国石油大学学报 (自然科学版) , 2007.2, 31 (1) :19.

[9]Collino F and Tsogka C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic hetero geneous media.Geophysics, 2001, 66 (1) :294~307.

[10]夏凡, 董良国, 马在田.三维弹性波数值模拟中的吸收边界条件[J].地球物理学报, 2004.1, 47 (1) :132.

高阶微分方程篇6

时域有限差分方法[finite-difference time-domain (FDTD) method]自1966年由Yee[1]提出来后, 一直是计算电磁学常用的方法之一, 并广泛应用于电磁散射、天线的分析与设计、雷达截口的计算等电子工业与国防工业。但是, 由于FDTD方法是条件稳定的, 即在二维情形下, 时间和空间步长分别为Δt, Δx, Δy;必须满足Courant-Friedrichs-Lewy (CFL) 稳定条件: cΔt≤[1/ (Δx) 2+1/ (Δy) 2]-1/2。其中c是在介质中的光速, ε, μ是介电常数和磁导率。为了解决上述问题, 由文献[2,3]提出了交替方向隐式时域有限差分方法 (ADI-FDTD) 方法, 并用Fourier方法证明了这种格式是无条件稳定的。近年来[4], 将分裂算子方法与FDTD方法结合提出新颖而且简单的分裂算子的时域有限差分方法 (S-FDTD) , 与二维ADI-FDTD相比都是无条件稳定的二阶格式, 但是S-FDTD格式采用算子分裂技术, 所以格式简单, 计算时间短, 并且在模拟一种散射问题时, S-FDTD比ADI-FDTD更精确。麦克斯韦方程的高阶时域有限差分方法已有较多的研究工作[4,5,6,7,8,9], 但是这些方法都没有使用分裂算子技术。

本文将高阶差分与S-FDTD方法相结合提出一种新的方法, 即在分裂后的方程基础上对空间采取四阶中心差分, 从而得到分裂 (splitting) 的高阶 (high-order) 的时域有限差分 (FDTD) 格式SHO-FDTDⅠ;并在此格式的基础上通过增加扰动项减小分裂误差给出它的修正格式SHO-FDTDⅡ, 给出了具体的计算步骤, 然后用Fourier方法证明了这两种格式是无条件稳定的, 并给出这两种格式的数值弥散关系式和数值弥散误差的计算。与条件稳定的高阶FDTD相比, 新方法是无条件稳定的;与高阶的ADI-FDTD相比, 新格式更简单, 易于编程和实现。

1 二维麦克斯韦方程的SHO-FDTD格式

考虑以下二维横向电波:

$\frac{\partial E_{x}}{\partial t} = \frac{1}{ε} \frac{\partial Η_{z}}{\partial y}$ (1.1)

$\frac{\partial E_{y}}{\partial t} = - \frac{1}{ε} \frac{\partial Η_{z}}{\partial x}$ (1.2)

$\frac{\partial Η_{z}}{\partial z} = \frac{1}{μ} (\frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x})$ (1.3)

其中 $\vec{E} = (E_{x} (x, y, t), E_{y} (x, y, t))$ 表示电场, $Η_{z} = Η_{z} (x, y, t)$ 表示磁场, ε和μ分别是介电常数和磁导率, Ω=[0, a]×[0, b], t ∈ (0, T]。其边界满足理想导体边界条件: $E_{x} (x, 0, t) = E_{x} (x, b, t) = E_{y} (0, y, t) = E_{y} (a, y, t) = 0$ 。初始条件: Ex0=Ex0 (x, y) , Ey0=Ey0 (x, y) , Hz0=Hz0 (x, y) 。

为了书写的方便, 我们仅讨论常系数的情形, 这里所用的方法容易推广到变系数的情形。对空间区域Ψ以及时间区域[0, T]采取如下剖分:

其中Δx和Δy分别是沿x轴方向和沿y轴方向的空间离散步长, Δt是时间步长, I, J, N为整数。

对于任意一个网格函数F (t, x, y) , 引入以下记号:

F $_{α, β}^{m}$ =F (mΔt, αΔx, βΔy) , $δ_{t} F_{α, β}^{m} = \frac{F_{α, β}^{m + \frac{1}{2}} - F_{α, β}^{m - \frac{1}{2}}}{Δ t},$

$δ_{x}^{2} F_{α, β}^{m} = \frac{1}{Δ x} (\frac{1}{24} F_{α - \frac{3}{2}, β}^{m} - \frac{9}{8} F_{α - \frac{1}{2}, β}^{m} + \frac{9}{8} F_{α + \frac{1}{2}, β}^{m} - \frac{1}{24} F_{α + \frac{3}{2}, β}^{m})$ ,

$δ_{y}^{2} F_{α, β}^{m} = \frac{1}{Δ y} (\frac{1}{24} F_{α, β - \frac{3}{2}}^{m} - \frac{9}{8} F_{α, β - \frac{1}{2}}^{m} + \frac{9}{8} F_{α, β + \frac{1}{2}}^{m} - \frac{1}{24} F_{α, β + \frac{3}{2}}^{m})$ 。

用E $_{v_{α, β}}^{m}$ (v=x, y) 表示电场Ev (tm, xα, yβ) 的近似值。H $_{z_{α, β}}^{m}$ 表示磁场Hz (tm, xα, yβ) 的近似值。将分裂算子的时域有限差分方法[9]与高阶差分方法[4]相结合, 提出如下格式:

第一步:

$\frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n})$ (1.4a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.4b)

第二步:

$\frac{E_{x_{i + \frac{1}{2}, j}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} + Η_{z_{i + \frac{1}{2}, j}}^{n})$ (1.5a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}}{Δ t} = \frac{1}{2 μ} δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.5b)

此格式称为分裂算子 (splitting) 的高阶 (high-order) 时域有限差分方法, 记为SHO-FDTDⅠ。其中, 边界条件为: $E_{x_{i + \frac{1}{2}, 0}}^{n} = E_{x_{i + \frac{1}{2}, j}}^{n} = E_{y_{0, j + \frac{1}{2}}}^{n} = E_{y_{i, j + \frac{1}{2}}}^{n} = 0$ 。初始值:

求解步骤为:先解第一步, 由式 (1.4b) 得:

$Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n} - \frac{Δ t}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.6)

将式 (1.6) 代入式 (1.4a) 整理得:

显然这是一个七对角的方程组, 系数矩阵元素都是常数而且在每一个时间层上都不变。因此可由一些求解线性方程组的方法如:超松弛迭代法, 共轭梯度法等, 求出 ${E_{y_{i, j + \frac{1}{2}}}^{n + 1}}$ 。然后, 将 $E_{y_{i, j + \frac{1}{2}}}^{n + 1}$ 代入式 (1.6) 显式求出 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。式 (1.5) 的计算与第一步相似, 但要用到第一步算出的 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。方程 (1.7) 中的第一个方程和最后一个方程的系数需要调整, 可根据文献[5]中的方法进行, 这里略去。

检查格式SHO-FDTDⅠ的截断误差发现关于时间它的精度不高。为了提高精度, 引入一个扰动项, 得到修正格式SHO-FDTDⅡ如下:

第一步:

第二步与SHO-FDTDⅠ中的第二步相同.

这种格式的初边值条件与式 (1.4) 、式 (1.5) 的初边值条件完全相同, 求解方法也与前一格式相同。为了了解格式SHO-FDTDⅡ的逼近精度, 由式 (1.5a) 、式 (1.5b) 、式 (1.8a) 、式 (1.8b) 消去中间项 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ , 可得到SHO-FDTDⅡ的等价格式:

$δ_{t} E_{x_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} = \frac{1}{2 ε} δ_{y}^{2} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{n}) (1.9 a)$

$\begin{array}{l} δ_{t} E_{y_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} = - \frac{1}{2 ε} δ_{x}^{2} (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n}) + \\ \frac{Δ t}{4 μ ε} δ_{x}^{2} δ_{y}^{2} (E_{x_{i, j + \frac{1}{2}}}^{n + 1} - E_{x_{i, j + \frac{1}{2}}}^{n}) (1.9 b) \end{array}$

$\begin{array}{l} δ_{t} Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = \frac{1}{2 μ} [δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) - \\ δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})] (1.10) \end{array}$

格式SHO-FDTDⅠ的等价格式与式 (1.9a) 、式 (1.9b) 、式 (1.10) 相似, 只要将式 (1.9) 最右端后一项中的“-”改为“+”。由此可以看出, SHO-FDTDⅡ的摄动误差 (二阶) 比SHO-FDTDⅠ的摄动误差 (二阶) 要高一阶, SHO-FDTDⅡ是关于时间2阶, 关于空间是4阶的式 (2.4) 格式, SHO-FDTDⅠ则是式 (1.4) 格式。

2 稳定性分析与数值弥散分析

用Fourier方法分析这两种格式的稳定性。假定格式的差分解具有下列形式。

$\vec{E}_{α, β}^{n} = \vec{E}_{0} ξ^{n} e^{- i (k_{x} α Δ x + k_{y} β Δ y)}$ ;

H $_{z_{α, β}}^{n}$ =Hz0ξne-i (kxαΔx+kyβΔy) 。

其中 $i = \sqrt{- 1}$ , $\vec{E}_{0} = (E_{x_{0}}, E_{y_{0}})^{Τ}$ 和Hz0是振幅, $\vec{k} = (k_{x}, k_{y})$ 是波矢量, $| k | = \sqrt{k_{x}^{2} + k_{y}^{2}}$ , ξ是增长因子。下面将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入式 (1.5a) 、式 (1.9) 、式 (1.10) , 化简得:

$(ξ - 1) E_{x_{0}} - i \frac{Δ t}{ε} (ξ + 1) b_{y} Η_{z_{0}} = 0$ (2.1)

$\begin{array}{l} (ξ - 1) E_{y_{0}} + \frac{(Δ t)^{2}}{ε μ} (ξ - 1) a_{x} b_{y} E_{x_{0}} + \\ i \frac{Δ t}{ε} (ξ + 1) a_{x} Η_{z_{0}} = 0 (2.2) \end{array}$

$\begin{array}{l} (ξ - 1) Η_{z_{0}} - i \frac{Δ t}{μ} (ξ + 1) b_{y} E_{x_{0}} + \\ i \frac{Δ t}{μ} (ξ + 1) a_{x} E_{y_{0}} = 0 (2.3) \end{array}$

其中 $a_{x} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{x} Δ x - \frac{9}{4} \sin \frac{1}{2} k_{x} Δ x)}{2 Δ x}, b_{y} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{y} Δ y - \frac{9}{4} \sin \frac{1}{2} k_{y} Δ y)}{2 Δ y}$ 。由于向量 (Ex0, Ey0, Hz0) 是非零向量, 所以关于Ex0, Ey0, Hz0的方程组的系数矩阵的行列式的值为零。整理得:

(ξ-1) (d0ξ2+2d1ξ+d0) 2=0 (2.4)

式 (2.4) 中

$d_{0} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2};$

$d_{1} = - 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}$ 。

解方程式 (2.4) 得: $ξ_{1} = 1, ξ_{2} = d_{0}^{- 1} (- d_{1} + i \sqrt{d_{0}^{2} - d_{1}^{2}}) 。 ξ_{3} = d_{0}^{- 1} (- d_{1} - i \sqrt{d_{0}^{2} - d_{1}^{2}})$ 。显然方程三个根的模都是1, 所以格式SHO-FDTDⅡ是无条件稳定的。对于格式SHO-FDTDⅠ, 同理可得:

c3ξ3+c2ξ2+c1ξ+c0=0 (2.5)

式 (2.5) 中

$\begin{array}{l} c_{3} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{2} = - 3 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}; \end{array}$

$\begin{array}{l} c_{1} = 3 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{0} = - 1 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2} 。 \end{array}$

方程 (2.5) 的根的表达式非常复杂 (为了简单, 这里省略) , 通过对根的表达式的研究发现 $| ξ | = 1 + Ο (Δ t)$ , 因此格式SHO-FDTDⅠ是耗散的无条件稳定的。为了直观的感受 $| ξ |$ , 下面给出它在不同情况下的变化曲线图。设kx=kcos φ。ky=ksin φ。其中k, φ。是向量k的圆柱坐标, 再设 $Ν_{λ} = \frac{λ}{h}, ω = c k, ω$ 是频率, λ是波长, Δx=Δy=h, 则Nλ是一个波长内的节点个数。

令 $S = \frac{c Δ t}{h}$ , 由CFL的定义可知

$c Δ t \sqrt{\frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}}} = \frac{c Δ t}{h} \sqrt{2} = \sqrt{2} S$ 。因此, 不妨称S为CFL数。由方程 (2.4) 或式 (2.5) 及其系数的表达式可知, 增长因子ξ是关于S, φ, Nλ的函数。例如, 可推出 $\sin \frac{3}{2} k_{x} Δ x = \sin (\frac{3 π}{Ν_{λ}} \cos φ)$ 。方程 (2.5) 的根的表达式比较复杂。下面我们用Matlab求出它的近似解, 并画出不同情况下它的根的变化曲线。

图1是在S=0.35, Nλ=10情况下给出了根 $| ξ |$ 随角度φ的变化曲线。上边的曲线是方程复数根的模, 下边的曲线是实数根的绝对值。很容易看出复数根的模比1大, 但是 $| ξ | = 1 + Ο (Δ t)$ , 其中, $Δ t = S \times \frac{1}{Ν_{λ}} = 0.035$ 。从图1中显然可以看出曲线位于: y=1±0.035两条直线之间, 这就表明SHO-FDTDⅠ格式是无条件稳定的。

图2和图3分别是在S=1.5, φ=35°和Nλ=10, φ=65°情况下, 方程 (2.5) 的根的模分别随着Nλ和S的变化情况, 从这两个图中同样可以看出SHO-FDTD格式是无条件稳定的。

2.2 数值弥散关系

假设麦克斯韦方程的差分解为: $\vec{E}_{α, β}^{n} =$ $\vec{E}_{0} e^{i (k_{x} α Δ x + k_{y} β Δ y - ω n Δ t)}$ , H $_{z_{α, β}}^{n}$ =Hz0ei (kxαΔx+kyβΔy-ωnΔt) , 其中, ω, kx, ky, $\vec{E}_{0}$ 和Hz0与2.1节的符号相同。将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入SH-FDTDⅡ的等价格式 (1.5a) 、式 (1.9) 、式 (1.10) , 类似于2.1节对增长因子ξ的推导, 可得:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2}] \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.6) \end{array}$

式 (2.6) 中 $c^{2} = \frac{1}{ε μ}$ , 式 (2.6) 就是SHO-FDTDⅡ的数值弥散关系式。

同理可得SHO-FDTDⅠ格式的数值弥散关系式:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2} \times \\ \cos (\frac{1}{2} ω Δ t) \sin^{- 1} (\frac{1}{2} ω Δ t)] \times \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.7) \end{array}$

由于 $\lim_{x \to 0} a_{x} = \lim_{x \to 0} a_{x}^{2} = \frac{k_{x}^{2}}{4}, \lim_{x \to 0} b_{y} = \lim_{x \to 0} b_{y}^{2} = \frac{k_{y}^{2}}{4}$ ;所以当Δx, Δy, Δt趋于0时, 数值弥散关系式 (2.6) 、式 (2.7) 趋向于理想的数值弥散关系: ω2=c2[k $_{x}^{2}$ +k $_{y}^{2}$ ]。此外, 比较式 (2.6) 、式 (2.7) 容易看出SHO-FDTDⅠ和SHO-FDTDⅡ的格式在时间上的精度分别是一阶的和二阶的, 并且SHO-FDTDⅡ比SHO-FDTDⅠ的数值弥散误差要小得多。

2.3 数值弥散误差

下面通过试验求出格式SHO-FDTDⅠ和格式SHO-FDTDⅡ的数值弥散误差, 并与理论结果进行对比分析。令ξ=eiωΔt是方程式 (2.1) 、式 (2.4) 的根。ω=ωR+iωI, 则ξ=e-ωIΔt[cos (ωRΔt) +isin (ωRΔt) ]从而, $\tan (ω_{R} Δ t) = \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ ; 其中Im (ξ) , Re (ξ) 分别表示ξ的虚部与实部。波的数值相速vp与光速c的比的数值弥散误差: $\frac{v_{p}}{c} = \frac{ω_{R}}{k c} = \frac{1}{c k Δ t} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)} = \frac{Ν_{λ}}{2 π S} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ , 其中S是CFL数, Nλ是一个波长内的节点数, k是波长值。

图4—图6给出了Vp/c在不同情况下的变化曲线。图4给出了SHO-FDTDⅠ与SHO-FDTDⅡ的Vp/c在S=0.35, Nλ=10的条件下随着角度φ的变化图, 从图4中可以看出SHO-FDTDⅡ格式的Vp/c更接近于1。图5给出了两种格式在S=2.4, φ=120°的条件下Vp/c随着Nλ的变化图, 从图5中可以看出SHO-FDTDⅡ格式的数值弥散误差小于SHO-FDTDⅠ的误差。图6是两种格式在φ=65°, Nλ=10的条件下Vp/c随着S的变化图, 从图6中可以看出随着S的增大, 数值弥散误差变得越来越大, 但是在所有的情况下, SHO-FDTDⅡ格式的数值弥散误差比SHO-FDTDⅠ格式的误差要小得多。

参考文献

[1] Yee K S.Numerical solution of initial boundary value problems invol-ving Maxwell’s equations in isotropic media.IEEE Trans Antennasand Propagation, 1966;14:302—307

[2] Zheng F, Chen Z, Zhang J.Toward the development of a three-dimen-sional unconditionally stable finite—difference time—domain method.IEEE Trans Microwave Theory Tech, 2000;48:1550—1558

[3] Namiki T.AnewFDTD algorithmbased on alternatingdirection implicitmethod.IEEE Trans Microwave Theory Tech, 1999;47:2003—2007

[4] Turkel E, Yefet A.Fourth order compact method for the Maxwell e-quations with discontinuous coefficients.Applied Numerical Mathe-matics, 2000;33:125—134

[5] Shang J S.Higher-order compact difference schemes for time depend-ent Maxwell’s equations.Journal of Computational Physics, 2003;153:312—333

[6] Xu L J, Yuan N C.高阶ADI-FDTD算法的数值色散分析.电子信息学报, 2005;27 (3) :1662—1665

[7] Xie Z Q, Zhang B, Chan C H.An explicit fourth-order staggered finitedifference time-domain method for Maxwell’s equations.Journal ofComputational and Applied Mathematics, 2002;147:75—98

[8] Manry C W, Broschat S L, Scneider J B.Higher-order FDTD meth-ods for large problems.J Applied Computational Electromagnetics So-ciety, 1995;10:17—29

【高阶微分方程】推荐阅读：

高阶P-Laplace方程边值问题的上下解方法09-23

高阶学习06-03

高阶能力07-29

高阶非线性06-03

高阶有限元10-11

培养学生高阶思维能力08-10

scala入门09 scala高阶函数10-13