差分方程

2024-10-29

差分方程（精选9篇）

差分方程篇1

对于常系数线性差分方程而言, 它的解法已经相当完善了.而对于非线性差分方程而言, 它的类型广, 形式多样, 所以并没有固定的解法.本文的主要思想就是通过各种变形将这类非线性差分方程转变为线性差分方程求解.

定理1 差分方程an=λan-1+ban-2的特征方程为x2-λx-b=0.若其特征根x1, x2为一对相异实根, 则差分方程的通解为an=c1x $_{1}^{n}$ +c2x $_{2}^{n}$ .

若x1=x2, an= (c1+c2n) x $_{1}^{n}$ .

若x1=δ+iω, x2=δ-iω, an=c1ρncosnθ+c2ρnsinnθ.

其中 $ρ = \sqrt{δ^{2} + ω^{2}} ‚ θ = \arctan \frac{ω}{δ} .$

定理2 差分方程 $a_{n} = λ + \frac{b}{a_{n - 1}}$ 的特征方程为x2-λx-b=0.若其特征根x1, x2为一对相异实根, 则差分方程的通解为 $a_{n} = \frac{c_{1} x_{1}^{n + 1} + c_{2} x_{2}^{n + 1}}{c_{1} x_{1}^{n} + c_{2} x_{2}^{n}} .$

若 $x_{1} = x_{2} ‚ a_{n} = \frac{[c_{1} + c_{2} (n + 1)] x_{1}^{n + 1}}{(c_{1} + c_{2} n) x_{1}^{n}} = \frac{c_{1} + c_{2} (n + 1)}{c_{1} + c_{2} n} x_{1} .$

$\begin{array}{l} 若 x_{1} = δ + i ω ‚ x_{2} = δ - i ω ‚ \\ a_{n} = \frac{c_{1} ρ^{n + 1} \cos (n + 1) θ + c_{2} ρ^{n + 1} \sin (n + 1) θ}{c_{1} ρ^{n} c o s n θ + c_{2} ρ^{n} s i n n θ} . \end{array}$

其中 $ρ = \sqrt{δ^{2} + ω^{2}} ‚ θ = \arctan \frac{ω}{δ} .$

证明令 $a_{n} = \frac{b_{n + 1}}{b_{n}}$ , 代入原差分方程, 得bn+1=λbn+bbn-1.应用常系数线性差分方程的解法即可得结论.

为简化计算, 设差分方程 $a_{n} = \frac{a_{n - 1}}{2} + \frac{k}{a_{n - 1}} (k > 0)$ 中 $k = \frac{1}{2} .$

令 $a_{n} = \frac{b_{n}}{c_{n}}$ 代入上式.

$\begin{array}{l} 令 {\begin{cases} b_{n} = b_{n - 1}^{2} + c_{n - 1}^{2} ‚ \\ c_{n} = b_{n - 1} c_{n - 1} ‚ \end{cases} \\ ∴ b_{n} + c_{n} = (b_{n - 1} + c_{n - 1})^{2} ‚ \\ ∴ b_{n} + c_{n} \geq 0 (n \geq 2) . \end{array}$

假设b1+c1≥0, ∴bn+cn≥0 (n≥1) .

应用线性差分方程的解法bn+cn= (b1+c1) 2n-1.

$\begin{array}{l} 同理 b_{n} - c_{n} = (b_{1} - c_{1})^{2^{n - 1}} ‚ \\ 2 b_{n} = (b_{1} + c_{1})^{2^{n - 1}} + (b_{1} - c_{1})^{2^{n - 1}} ‚ \\ 2 c_{n} = (b_{1} + c_{1})^{2^{n - 1}} - (b_{1} - c_{1})^{2^{n - 1}} ‚ \\ a_{n} = \frac{2 b_{n}}{2 c_{n}} = \frac{(b_{1} + c_{1})^{2^{n - 1}} + (b_{1} - c_{1})^{2^{n - 1}}}{(b_{1} + c_{1})^{2^{n - 1}} - (b_{1} - c_{1})^{2^{n - 1}}} . \end{array}$

$\begin{array}{l} 令 c_{1} = 1 ‚ b_{1} = a_{1} ‚ \\ a_{n} = \frac{(a_{1} + 1)^{2^{n - 1}} + (a_{1} - 1)^{2^{n - 1}}}{(a_{1} + 1)^{2^{n - 1}} - (a_{1} - 1)^{2^{n - 1}}} . \end{array}$

对于任意常数k, 差分方程 $a_{n} = \frac{a_{n - 1}}{2} + \frac{k}{a_{n - 1}} (k > 0) .$

$\begin{array}{l} 令 a_{n} = \sqrt{2 k} b_{n} ‚ a_{n - 1} = \sqrt{2 k} b_{n - 1} ‚ \\ ∴ b_{n} = \frac{b_{n - 1}}{2} + \frac{1}{2 b_{n - 1}} ‚ \\ b_{n} = \frac{(b_{1} + 1)^{2^{n - 1}} + (b_{1} - 1)^{2^{n - 1}}}{(b_{1} + 1)^{2^{n - 1}} - (b_{1} - 1)^{2^{n - 1}}} \\ = \frac{(\frac{a_{1}}{\sqrt{2 k}} + 1)^{2^{n - 1}} + (\frac{a_{1}}{\sqrt{2 k}} - 1)^{2^{n - 1}}}{(\frac{a_{1}}{\sqrt{2 k}} + 1)^{2^{n - 1}} - (\frac{a_{1}}{\sqrt{2 k}} - 1)^{2^{n - 1}}} ‚ \\ ∴ a_{n} = \sqrt{2 k} b_{n} \\ = \sqrt{2 k} \cdot \frac{(\frac{a_{1}}{\sqrt{2 k}} + 1)^{2^{n - 1}} + (\frac{a_{1}}{\sqrt{2 k}} - 1)^{2^{n - 1}}}{(\frac{a_{1}}{\sqrt{2 k}} + 1)^{2^{n - 1}} - (\frac{a_{1}}{\sqrt{2 k}} - 1)^{2^{n - 1}}} . \end{array}$

定理3 差分方程 $a_{n} = \frac{a_{n - 1}}{2} + \frac{k}{a_{n - 1}} (k > 0)$ 的通解为:

参考文献

徐利治.计算组合数学[M].上海:上海科学技术出版社, 1983:26-28.

差分方程篇2

从解的渐近状态着手,将所述方程的所有非平凡解分成互不相交的`四类,应用分类讨论方法和分析方法,讨论了一类广泛的二阶非线性差分方程解的振动性质,建立了两个新的振动性定理,推广并改进了已有文献中的相关结果.

作者：张全信燕居让 ZHANG Quan-xin YAN Ju-rang 作者单位：张全信,ZHANG Quan-xin(滨州学院数学与信息科学系,山东滨州,256603)

燕居让,YAN Ju-rang(山西大学数学科学学院,太原,030006)

差分方程篇3

自动聚焦完成的优劣,是成像系统能否获得清晰图像的最关键一步。自动聚焦作为核心技术之一,是走向数字产品品牌化、自主化、加强市场竞争力的关键技术之一; 也是人类突破人眼自身生理局限,实现数字产品的智能化应用的关键技术之一。它的应用已经超越了传统的数码相机、元器件检验、医学成像、空间探测、显微成像,而且延伸到了安全防卫、远程医疗、机器视觉等领域[1,2]。动聚焦算法的改进是在分析传统峰值搜索算法和聚焦曲线特点的基础上进行的。将聚焦曲线附近的区域取名为聚焦区,将远离焦点的区域取名为散焦区。通过分析得知,聚焦曲线在散焦区和聚焦区分别呈现出各自的特点。本文采用TI的DM368 处理器的硬件实现的聚焦评价函数,在散焦区采用基于差分方程预测模型的峰值搜索算法预测峰值搜索的方向,在聚焦区采用单调性判断法搜索峰值对应的聚焦位置,实现了自动聚焦系统。通过实验表明,本文改进的算法有效地解决了视频质量反复变化、峰值搜索时间长和散焦区局部极值干扰的问题。

1 传统峰值搜索算法

峰值搜索算法是一个一维数列的极值最优化问题,是对深度优先搜索的一种改进。在自动聚焦算法中,峰值搜索算法的关键是对聚焦方向的控制,确保快速而准确地搜索到聚焦曲线的峰值,即聚焦位置。在自动聚焦领域,最常用到的峰值搜索算法有: 单调性判断法、斐波那契搜索算法、曲线拟合算法、盲人爬山算法等。

传统的单调性判断法是通过两个值比较来判断此时所处聚焦曲线的位置,从而寻找下一步的搜索方向。此方法适用于理想的没有局部最大值的曲线,抗干扰能力差,容易导致错误聚焦[3]。

盲人爬山算法是推动电机从负轴向正轴进行一次遍历,遍历结束后返回至聚焦评价函数输出的最大值对应的电机位置[4,5]。盲人爬山算法可以有效地消除局部最大值的干扰,在静态摄像机中可以达到聚焦效果。但是,在实际应用中,一次遍历需要消耗大量的聚焦时间,无法满足高清网络摄像机聚焦实时性的要求。

斐波那契搜索算法来源于斐波那契数列的概念。斐波那契搜索算法是通过缩小搜索区间的范围来搜索单峰曲线极点的[6]。但是,斐波那契搜索算法只是一种理论计算上的最优化单峰搜索算法,实际的应用中存在诸多的限制条件和缺点。

本文提出的峰值搜索算法是基于差分预测模型的峰值搜索算法与传统的单调性判断法相结合的方法。通过充分地分析聚焦曲线的特点,在聚焦区,聚焦值变化量大,聚焦曲线陡峭,单调性十分明显,不会出现局部最大值,采用传统的单调性判断法;在散焦区,聚焦值变化量小,聚焦曲线波动较大,局部极值现象普遍存在,采用差分方程预测模型的峰值搜索算法。

2 聚焦曲线特点分析

聚焦值和聚焦位置的关系称作聚焦值曲线。横轴表示聚焦位置,纵轴表示聚焦值,如图1 所示。

聚焦值曲线最大值对应的聚焦位置为焦点。远离焦点的区域为散焦区,焦点附近的区域为聚焦区。散焦区与聚焦区的判定: 判断聚焦值的变化量是否大于设定值,大于则认为进入聚焦区; 反之,进入散焦区。

通过分析可知聚焦曲线有如下的特点:

1) 在聚焦区,聚焦值曲线非常陡峭; 在散焦区,聚焦值曲线较为平坦。

2) 在散焦区,聚焦值曲线有明显波动,局部最大现象较为明显; 在聚焦区,聚焦值曲线波动较小,但是跨度较小。

根据图1 特点分析,设计峰值搜索算法: 在散焦区,采用基于差分方程预测模型的峰值搜索算法进行搜索方向的确定; 在聚焦区,采用单调性判断法进行峰值搜索。

3 差分方程预测模型

灰色预测模型是通过少量的、不完全的信息,建立数学模型并做出预测预测的一种预测方法[7,8,9]。灰色预测模型是在假定生成序列满足微分方程的基础上进行的,而差分方程预测模型是在假定生成序列满足差分方程的基础上进行预测的。通过实验表明,使用差分方程的预测,预测曲线拐点出现的位置比灰色预测模型更加准确。

1)差分方程预测模型原理框图

GM(1,1)原理如图2所示。

此处,x( 0)代表原始的数据序列。AGO表示累加生成,即通过数列间各时刻数据的依个累加以得到新的数据。在建立差分方程预测模型之前,需先对原始时间序列进行数据处理,经过数据处理后的时间序列即称为生成列。经过累加生成弱化了原始时间序列的随机性,强化规律性,对原始数据的生成就是企图从杂乱无章的现象中去发现内在规律。DDE( 2,1) 表示差分方程模型。IAGO表示累减生成,是累加生成的逆运算。累减生成,即对数列求相邻两数据的差,累减生成是累加生成的逆运算,常简记为IAGO,累减生成可将累加生成还原为非生成数列。

2)差分方程模型的建立详细步骤

①获取n个原始数据序列

式中,x( 0)( n) 表示第n个原始数据。

② 运用AGO运算,得到累加生成序列

其中。

③ 建立差分方程预测模型

其中,a和b为差分方程待确定的系数,t为整数。

为了求得差分方程的系数,使用最小二乘法。差分方程可以表述为,方程可以表述为:

使p = 1,2,…,n - 1 ,可以得到方程:

可以使用线性最小二乘法,得到估值参数:

令x( 1)( p) = rp,求得其对应的特征方程:

解之可得:

3) 差分方程求解

根据特征根的不同情况,求解方程式( 3) 的解,求解过程如下:

① 若方程式( 3) 有两个互不相同的实根r1和r2,则其解为:

考虑到初始情况下p = 1,p = 2 时,可得到方程式( 9) 和方程式( 10) :

求解方程组式( 9) 和式( 10) ,可以得到常数C1和C2如下:

② 若方程式( 3) 有两个相同的实根,则其解为:

考虑到初始情况下p = 1,p = 2 时,可得到方程式( 12) 和方程式( 13) :

求解方程组式( 12) 和式( 13) ,可以得到常数C1和C2如下:

③ 若方程式( 3) 的根为复数,则其解为:

其中:

考虑到初始情况下p = 1,p = 2 时,可得到方程式( 15) 和方程式( 16) :

求解方程组式( 15) 和式( 16) ,可以得到常数C1和C2如下:

4 聚焦步长确定

若以恒定速度进行聚焦,则存在如下问题。若设定步长太大( 即控制镜头移动的时间太长) ,则会出现聚焦错误或找不到焦点。因为聚焦区的距离很短。若设定的步长太小,则会出现聚焦时间太长的问题[10]。分析聚焦曲线可知,步长的大小与聚焦值变化量应该成反比关系。

在散焦区采用大步长进行搜索,在聚焦区采用小步长进行搜索。在聚焦过程中,为了保证快速而准确地找到焦点对应的位置,需要实时的调整搜索步长。

在散焦区,使用较大步长,从而保证预测方向的准确性。搜索策略如图3 所示。

步长调整策略如下:

步骤一

设定初始步长为1,连续采样聚焦值三次,判断当前所在区域。若为散焦区,跳至步骤二; 若为聚焦区,使用爬山算法进行峰值搜索。

步骤二

增加步长为2,采集新的聚焦值,判断聚焦值的的变化量,确定当前所在区域。若为散焦区,使用差分方程预测模型预测搜索方向,跳至步骤三。若为聚焦区,使用最小步长的单调性判断法进行峰值搜索。

步骤三

增加步长为4,进行聚焦值的采集,判断当前区域并进行搜索方向的预测。若为聚焦区,改用最小步长进行峰值搜索。若为散焦区,继续采集聚焦值,进行下一步的判断。

在散焦区搜索峰值过程中,搜索步长不宜过长,否则有可能直接跨过聚焦区[11,12],也可能造成图像质量出现明显的反复变化。在聚焦区,使用较小步长进行搜索。这样既快速地找到焦点对应的位置,又能避免图像质量的明显变化。

5 实验结果分析

5. 1 硬件平台设计

本实验的硬件实验平台设计如图4所示。

信号转换模块:表示一个串口转485的信号转换电路。球机表示一台带有红外智能球型云台的高清网络摄像机。实验板代表基于TMS320DM368处理器的高清网络摄像机的实验板,主要功能是利用处理器的AF引擎作为聚焦评价函数的输出。

① 表示PC端应用软件的输出信号,用来控制红外智能球型云台。② 表示信号转化模块的输出信号,用来控制红外智能型云台。③ 表示红外智能球型云台视频输出的视频信号。④表示PC向基于TMS320DM368 处理器的高清网络摄像机的实验板发出的控制命令信号。⑤ 表示基于TMS320DM368 处理器的高清网络摄像机的实验板向PC端发送的信息。

5. 2 应用软件设计

应用软件: 使用Microsoft Visual C#语言编写,采集网口的数据进行运算,根据运算结果发送聚焦命令,最终通过RS232 接口控制红外智能球型云台。软件界面如图5 所示。

差分方程预测模型: 使用Matlab R2011a实现,然后通过Matlab Deploy Tool生成COM组建,供应用应用程序调用。聚焦控制命令: 基于派尔高-D协议。云台左转向和云台右转向控件,用来控制云台的转动方向,主要用来变换场景。聚焦速度慢、中和快控件,分别用来设置聚焦的速度。停止取数和接受数据分别表示将web页面取来的数据保存至txt文档中,供后台的分析和应用。聚焦近和聚焦远表示发送聚焦命令,手动控制聚焦的方向。自动聚焦控件用来触发自动聚焦。

5. 3 散焦区实验结果分析

在聚焦调节过程中,采用GM( 1,1) 模型的搜索算法可以准确地确定聚焦的方向,最终能快速地找到聚焦区。

聚焦步长为4,聚焦位置为{ 1,5,9,13,17,21,25,29} ,采集到的聚焦值为聚焦值分别为{ 24 810 403,24 811 703,25 155 444,25 227 081,25 281 798,25 669 023,25 668 664,26 440 896} 。散焦区聚焦采样如图6“* ”所示。

选取{ 24 810 403,24 810 703,25 155 444,25 227 081} 为初始值,初始方向为远,进行GM( 1.1) 模型试验。实验结果如表1 所示。

如表1 所示,在远方向的聚焦过程中,采用步长为4 的搜索,预测的聚焦方向都为远,可以保证准确地搜索到聚焦区。由后验差比值C和小误差概率分析可知,预测精度等级都达到了优秀级别。使用传统的爬山算法,则会受到25 669 023 点处局部最大值的影响。

通过较大量实验证明,GM( 1,1) 模型进行的方向预测的准确率可达90% 以上,而极限预测的时间不超过19. 1 ms,小于相邻两帧的间隔时间33. 3 ms,从而保证了聚焦的速度。

5. 4 聚焦区实验结果分析

在聚焦区,聚焦曲线非常地陡峭,在聚焦位置为{ 29,30,31,32,33,34 } 处,对应的聚焦值分别为{ 26 440 896,27 540 896,28 870 830,30 824 367,35 107 979,33 124 367} 。聚焦区聚焦采样如图5“o”所示。

在聚焦区,聚焦值的单调性很明显。使用小步长的爬山算法,能快速地找到峰值,从而聚焦成功。同时,可以避免图像质量的反复变化。

此外,聚焦过程中计算量很小。在此过程中容易受场景变化的影响,使用退回检测最大值的方法有效地消除了干扰。

6 结语

通过大量实验,TMS320DM368 的AF引擎计算出每帧聚焦值的平均时间为0. 95 ms,而相同实验环境下,取1000 帧视频,灰度方差函数平均需要81 ms,梯度函数平均需要101 ms。随着百万像素摄像机的普及,采用硬件实现的AF引擎将是不可扭转的趋势。

聚焦评价函数与峰值搜索算法平均总耗时为19. 2 ms,最大耗时为20. 2 ms,小于相邻两帧的时间33. 3 ms。满足聚焦实时性的要求。GM( 1,1) 模型很好解决了爬山算法局部最大值干扰的问题。提高了聚焦的准确性。

实验表明,本文改进的自动聚焦算法具有更好的准确性和实时性。

参考文献

[1]Zhang H,Zhang Z,Yang H,et al.Real-time auto-focus system design based on climbing algorithm and its FPGA implementation[C]//Computational Intelligence and Security(CIS),2012 Eighth International Conference on.IEEE,2012:332-335.

[2]王彦芳.自动聚焦系统中评价函数性能与动态区域选取的研究[D].山东:山东大学,2012:8-10.

[3]黄伟琼,游林儒,刘少君.基于改进的灰度对比度函数的自动对焦方法[J].计算机应用,2011,31(11):3008-3014.

[4]胡凤萍,常义林,马彦卓,等.视频自动聚焦的实现研究[J].光子学报,2010,39(10):1901-1906.

[5]Wang D,Yin S,Chen C,et al.Application of GM(1,1)Model on Predicating the Outpatient Amount[M]//Advances in Multimedia,Software Engineering and Computing Vol.2.Springer Berlin Heidelberg,2012:71-76.

[6]Gamadia M,Rahman M,Kehtarnavaz N.Performance Metrics for Passive Auto-Focus Search Algorithms in Digital and Smart-Phone Cameras[J].Journal of Imaging Science and Technology,2011,55(1):69-70.

[7]Dongchen Tsai,Homer H Chen.Smooth control of continuous autofocus[C]//Orlando,FL:19th IEEE International Conference on Image Processing(ICIP),2012:373-376.

[8]Quoc Kien Vuong,Jeongwon Lee.Initial direction and speed decision system for auto focus based on blur detection[C]//Las Vegas,NV:IEEE International Conference on Consumer Electronics(ICCE),2013:222-223.

[9]方建卫.基于离散差分方程模型人口结构预测[J].华章,2012(1):348.

[10]Tsai,Dongchen,Homer H.Autofocus method:America,8254774[P].2012-08-28.

[11]Yata,Kunio.Autofocus system:America,8363108[P].2013-01-29.

差分方程篇4

解KDV方程的一个二阶三层差分格式

本文讨论KDV方程定界问题的数值解法.对KDV方程构造了一个二阶三层的差分格式,并对非线性项进行了线性化,使格式的近似解更精确.通过严格的.误差估计证明了非线性稳定性.数值实验结果表明了理论证明的正确性和格式的有效性.该格式是可行的,有推广价值.

作者：李家永阿不都热西提・阿不都外力 LI Jia-yong Abudurexiti・Abuduwaili 作者单位：新疆大学数学与系统科学学院,乌鲁木齐,830046刊名：工程数学学报 ISTIC PKU英文刊名：CHINESE JOURNAL OF ENGINEERING MATHEMATICS年，卷(期)：24(6)分类号：O241.82关键词：KDV方程非线性项非线性稳定性误差估计

一类中立型差分方程解的振动性篇5

关键词：高阶时滞差分方程,振动,最终正解

众所周知, 差分方程理论在现代科技中是强有力的数学工具之一, 特别是生物种群动力学等前沿学科的迅速发展, 更促进了差分方程理论的再发展.尤其在离散系统方面, 描述突变现象等要比连续系统更高一筹.例如, 简单的模型对连续型其解是单调的, 而对差分方程其解就显得较复杂, 即随着参数的取值可能产生混沌和分支, Logistic模型就是如此.文献[3]研究差分方程的一般理论和工作, 文[5~6]研究了高阶差分方程解的振动性和渐近性.由于问题有较强的应用性, 相应的工作量和难度较大, 引起学者们的极大重视.

1. 振动性问题

考虑中立型差分方程:

其中d是奇数, τ>0, qi (n) >0 (i=1, 2, 3…) , σi是整数 (i=1, 2, 3…) .

向前差分Δ定义为Δxk=xk+1-xk且Δm=Δ (Δm-1) , (1.1) 的特殊情况已在泛函微分方程的数值分析中出现, 最近, 已有不少工作研究该方程的正解的存在的充要条件, 但对于该方程的解的振动性, 只是做了零星的工作.

本文将证明 (1.1) 的每个解都振动的充要条件是方程

下面, 若没有特别申明, 一个差分不等式是指对一切充分大的整数成立.

2. 方程 (1.1) 与方程 (1.2) 振动的等价性

设f (x) 和以下的f (x) 均是非减函数, 且xf (x) ≥0, f (x) =f (-x) , 则有以下的几条:

引理2.1方程的每个解振动的充要条件是差分不等式没有最终正解.

引理2.2方程 (m为偶数) 振动的充要条件是不等式Δmyn-1+pnf (yn) >0没有正解

证明参见文献[5].

证明不失一般性, 以d=3的情况给出证明.

(充分性) 设不然, 让{xn}为方程的最终正解, 设yn=xn-xn-τ, 则Δ3yn≤0可以证明, 最终地有yn>0, 因此ynΔ3yn≤0.由文献[1]定理1.7.11, 有两种情况:

(1) yn>0, Δyn<0, Δ2yn>0;

(2) yn>0, Δyn>0, Δ2yn>0.

首先考虑{xn}是 (1) 种情况的类型.

在这种情况下, , 设N是充分大的整数, 使当n≥N-τ时, 有xn≥0, yn≥0, Δyn≤0, Δ2yn≤0.

由于可令x=f (x) , 因此满足引理2.2, 方程有一正解, 这是一个矛盾.

当{xn}是第二种情况时, 可类推证明.

现就情形〈1〉来讨论

, 若l有限, 则k为常数, 于是存在正常数N, M, 当n≥N-2-min{σi} (i=1, 2, 3, …) 时有yn>M>0和Δyn

很明显Zn-Zn-1=Yn, 当n≥N时有Zn-Zn-1=Δyn-1

当n=N时有Zn=Δyn-1+Zn-1≤Δyn-1;

由引理2.1说明 (1.1) 有最终正解, 这是矛盾的, 利用文献[4]的方法容易证明 (2) 的情形.定理1证毕.

3. 应用举例

考虑如下中立型差分方程解的振动性:

其中α, β∈ (1, 2) 直接利用定理1和推论2就可以证明该方程的每一个解振动.

参考文献

[1]BRAYTON R K, WILLOUGHBY R A.On the numerical integration of a symmetric system ofdifference differential equation ofneutral type[J].JMathAnalApp, l, 1998 (18) :182-189.

[2]ZHANG Guang, CHENG Sui-sun.Oscillation criteria for a neutral difference equation delay[J].ApplMath Lett, 1995, 8 (3) :13-17.

[3]KELLEYW G, PETERSON A C.Difference equations:an introduction with applications[M].New York:Academic Press, 1991:1-288.

[4]AGARWALR P.Difference equation and inequalities[M].NewYork:MarcelDekker, 1990.

[5]张炳根, 杨博.非线性高阶差分方程的振动性[J].数学年刊, 1999, 20A (1) :71-80.

差分方程篇6

关键词：差分方程,平衡点,全局渐近稳定

1引言

文献[1]研究了高阶差分方程

xn+1=f (xn-s, xn-t) , n=0, 1, 2, …,

s>t, s, t∈{0, 1, 2, …}, (1)

其初始条件为x-s, x-s+1, …, x0∈ (0, +∞) , f满足:f (u, v) 关于u是递增的, 关于v是递降的.证明了这个结论:若f在I×I→I (令I=[a, b]) 上连续且满足:

(ⅰ) f (u, v) 在u上递增, 在v上递减;

(ⅱ) 当 (m, M) ⊂[a, b], 方程组有唯一解m=M.

则方程 (1) 有唯一平衡点 $\bar{x}$ , 且方程 (1) 的点都收敛于 $\bar{x} .$

受如上问题的启发, 本文考虑了定理中函数f的另一种情形, 并得出了有关结论.

2主要结论

定理令I=[a, b], 假设f在I×I→I上连续且满足:

(ⅰ) f (u, v) 在u上不增, 在v上不减;

(ⅱ) 当 (m, M) ∈[a, b], 方程组有唯一解m=M.

则方程 (1) 有唯一平衡点 $\bar{x}$ , 且方程 (1) 的点都收敛于 $\bar{x} .$

证明考察如下差分系统

${\begin{cases} m_{n} = f (Μ_{n - 1} ‚ m_{n - 1}) ‚ \\ Μ_{n} = f (m_{n - 1} ‚ Μ_{n - 1}) ‚ \end{cases}$

其中m0=a, M0=b.

由函数的单调性知

m1=f (M0, m0) ≥m0,

M1=f (m0, M0) ≤M0,

由数学归纳法进一步有:n=1, 2, …

a≤mn-1≤mn≤b,

a≤Mn≤Mn-1≤b.

由单调有界原理知数列{mn}{Mn}都收敛, 设 $\lim_{n \to \infty} m_{n} = m ‚ \lim_{n \to \infty} Μ_{n} = Μ$ , 则 (m, M) ⊂[a, b].由f的连续性知

因此 $m = Μ = \bar{x}$ , 则 $\bar{x} = f (\bar{x} ‚ \bar{x}) .$

即差分方程 (1) 有一个平衡点.

下证对任意自然数n有下面的式子成立

${\begin{cases} m_{n + 1} \leq x_{n (s + 1) + 1} \leq Μ_{n + 1} ‚ \\ m_{n + 1} \leq x_{n (s + 1) + 2} \leq Μ_{n + 1} ‚ \\ \dots \\ m_{n + 1} \leq x_{(n + 1) (s + 1)} \leq Μ_{n + 1} ‚ \end{cases} (2)$

利用数学归纳法证明.当k=0时, 由于

m1=f (M0, m0) =f (b, a)

≤x1=f (x-s, x-t)

≤f (a, b) =f (m0, M0) =M1,

m1=f (M0, m0) =f (b, a)

≤x2=f (x1-s, x1-t)

≤f (a, b) =f (m0, M0) =M1,

…

m1=f (M0, m0) =f (b, a)

≤xs+1=f (x0, xs-t)

≤f (a, b) =f (m0, M0) =M1.

即k=0时成立.

假设n=l时结论成立, 即就是

ml+1≤xl (s+1) +1≤ML+1,

ml+1≤xl (s+1) +2≤ML+1,

…

ml+1≤x (l+1) (s+1) ≤ML+1.

那么当n=l+1时, 由函数f的单调性以及数列{mn}{Mn}的单调性, 有

ml+2=f (Ml+1, ml+1) ≤x (l+1) (s+1) +1

=f (xl (s+1) +1, xl (s+1) +1+s-t)

≤f (ml+1, Ml+1) =Ml+2,

ml+2=f (Ml+1, ml+1) ≤x (l+1) (s+1) +1

=f (xl (s+1) +2, xl (s+1) +2+s-t)

≤f (ml+1, Ml+1) =Ml+2,

…

ml+2=f (Ml+1, ml+1) ≤x (l+2) (s+1)

=f (x (l+1) (s+1) +2, x (l+1) (s+1) +s-t)

≤f (ml+1, Ml+1) =Ml+2.

即当n=l+1时, 结论成立.

由 (2) 式及数列{mn}{Mn}的性质得 $\lim_{n \to \infty} x_{n} = \bar{x}$ .证毕!

参考文献

[1]樊永红, 王琳琳.一类高阶差分方程的全局渐近稳定性[J].兰州大学学报 (自然科学版) , 2008, (S1) :174-175.

[2]FEUER J, JANOWSKIEJ, LADAS G.Lyness-type equation in the third quadrant[J].Nonlin-ear Anal, 1997, 30 (2) :1183-1189.

[3]DeVault R, kosmala W, Ladas Gand Schultz SW.Global behavior ofyn+1= (p+yn-k) / (qyn+yn-k) [J].Nonliear Analysis TMA, 2001, 47:4743-4751.

差分方程篇7

关键词：时滞差分方程,特征方程,渐近稳定

1 引言

目前很多学者研究了有理型差分方程

一些特殊情形的稳定性,上述方程没有考虑时滞因素,在这篇文章中我们考虑上述方程的其中一类附加时滞项的情况,本文研究有理型时滞差分方程

的稳定性,得到了上述方程零解渐近稳定的充分条件,其中A≠0,C≠0,β4B4是常数,k>1正整数。

定理1.1设A≠0,C≠0,k>1正整数,则方程(1.1)零解渐近稳定的一个充分条件是下列条件之一成立:

设

准是方程

2 相关引理

考虑时滞差分方程

其中a,b为常数,k正整数,且ab≠0。

(2.0)的特征方程为

其中ab≠0。令,则方程(2.1)变为

其中,k=1,2,3或4(mod5)且c≠0。

记

引理2.1(i)k是奇数

(a)c>0时,方程(2.2)没有负根,特别地,0<c<β时,方程(2.2)有两个正实根,分别记为μ1(c)和μ2(c),它们分别属于(0,α)和(α,1),μ1(c)单调增加,μ2(c)单调减少;c=β时,α是方程(2.2)的一个二重根;c>β时,没有实根。

(b)c<0时,方程(2.2)有一个正根和一个负根,分别记为μ3(c)和μ4(c),它们分别属于(1,+∞)和(-∞,0),μ3(c)单调减少,μ4(c)单调增加。

(ii)k是偶数

(a)c>0时,方程(2.2)在(-∞,0)上有一个负根,记为μ5(c),其单调减少。特别地,0<c<β时,方程(2.2)分别在(0,α)和(α,1)上各有一个正根μ6(c)和μ7(c),且μ6(c)单调增加,μ7(c)单调减少;c=β时,α是方程(2.2)的二重根;c>β时,方程(2.2)没有实根。

(b)c<0时,方程(2.2)在(1,+∞)有一个正根记为μ8(c),其单调减少,且没有负根。

证明:记F(μ)=μk+5-μk+c,则,且F(0)=F(1)=c。下面我们只考虑k是奇数的情形,k是偶数的情形同理可证。c>0时,F(0)=c>0,F(-∞)=∞,且在(-∞,0)上F′(μ)<0,于是,方程(2.2)没有负根。特别地,0<c<β时,F(α)=c-β<0,F(0)=F(1)=c>0,且F′(μ)<0,μ∈(0,α)时,F′(μ)>0,μ∈(α,1)时,故方程(2.2)有两个正实根且分别位于(0,α)和(α,1)之间,记为μ1(c)和μ2(c);c=β时,F′(α)=0,而F′(α)≠0,故知α是一个二重根;c>β时,函数F(μ)的极小值F(α)=c-β>0,故方程(2.2)没有实根。至此得到(i)(a)。c<0时F(1)=c<0,F(∞)>0,且F′(μ)<0,故知当μ∈(1,+∞)时,方程(2.2)有唯一正根。又因为F(0)=c<0,F(-∞)>0,且在(-∞,0)上F′(μ)<0,故知当μ∈(-∞,0)时,方程(2.2)有唯一负根。这就得到(i)(b)。证毕。

考虑区间(0,π)的如下子区间,置

引理2.2设k是正整数,则

引理2.3当c增加时,方程(2.2)的根的绝对值也增加,除了下面几种情况:

(a)当k是奇数且0<c<β时,方程(2.2)较大的正实根随着c的增加而减小。

(b)当k是偶数且0<c<β时,方程(2.2)较大的正实根随着c的增加而减小。

(c)c>0,k≡1(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ1是方程G(θ)=0在上的根。

(d)c>0,k≡2(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ2是方程G(θ)=0在上的根。

(e)c>0,k≡3(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ3是方程G(θ)=0在π上的根。

(f)c>0,k≡4(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ4是方程G(θ)=0在,25π上的根。

(g)c<0,k≡1(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ5是方程G(θ)=0在上的根。

(h)c<0,k≡2(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ6是方程G(θ)=0在上的根。

(i)c<0,k≡3(mod5),则方程(2.2)辐角位于的复根的模随着c的增加而减小,这里θ7是方程G(θ)=0在上的根。

(j)c<0,k≡4(mod5),则方程(2.2)辐角位于的复根的模随着的增加而减小,这里θ8是方程上的根。

引理2.8对于线性系统来说,一致渐近稳定等价于指数稳定。

3 定理1.1的证明

容易知道

则方程(1.1)变为

其中

对(3.1)作自变量平移,则(3.1)可表成

其中,b=-C。我们知道方程(1.1)式与方程(3.3)式关于稳定性等价,又方程(2.0)为(3.3)的首次近似线性方程,于是,由(3.2)式,引理及定理,即得定理。证毕。

参考文献

[1]S.A.Kuruklis.The Asymptotic Stability of xn+1-axn+bxn-k=0,[J].Math.Anal.Appl.,1994,188:719-731.

[2]Hongshan Ren,Stability Analysis of Second Order Delay DifferenceEquations,[J].Funkcialaj Ekvacioj,2007,50:405-419.

[3]Yangtao Wang and Hongshan Ren,Asymptotic Stability of Third DelayDifferenceEquationswithTwoConstantCoefficients[J].Aduanceon Biomathematics,Volume II,Proceeding of the 6th Conference of Biomathematics,Liverpool,UK:World Academic press,2008:738-739.

[4]文斌,任洪善.一类四阶时滞差分方程的渐近稳定性研究,黑龙江大学自然科学学报,2010,27(3):323-331.

[5]F.M.Dannan,The Asymptotic Stability of xn+k+axn+bxn-l=0,[J].Journal ofDifference Equations and Applications,2004,10(6):589-599.

差分方程篇8

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.

差分方程篇9

时域有限差分方法[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

【差分方程】推荐阅读：

背景差分06-07

差分能量07-12

差分模型07-23

差分分析07-30

有限差分08-15

差分相关09-02

差分解调09-22

差分故障10-25

差分进化策略07-05

差分跳频09-07

>> 查看更多相关文档