基于波动方程有限差分算法的接收函数正演与偏移
王红落,常 旭,陈传仁.基于波动方程有限差分算法的接收函数正演与偏移.地球物理学报,2005,48(2):415~422Wang H L,Chang X,Chen C R.Receiver function forward modeli ng and migration based on wave -equation finite difference method.Chinese J .Geop h ys .(in Chinese),2005,48(2):415~422
基于波动方程有限差分算法的接收函数正演与偏移
王红落1
,常 旭1
,陈传仁
2
1中国科学院地质与地球物理研究所,北京 1000292长江大学地球物理与石油资源学院,荆州 434023
摘 要 针对接收函数正演与偏移,本文采用波动方程有限差分算法.借鉴成熟的勘探地震学方法,引入等效速度概念,建立接收函数转换波与地震勘探反射波的等效走时方程,实现了基于波动方程有限差分算法的接收函数正演与偏移.数值计算表明,波动方程有限差分叠后偏移方法可以对点绕射和穹隆构造模型实现高精度成像.本文利用数值计算讨论了波动方程有限差分叠后偏移与Kirchhoff 叠后偏移对于接收函数偏移的适用性,还对偏移过程中速度模型的误差进行了分析.
关键词 有限差分正演,接收函数阵列,共转换点(CCP)叠加,有限差分偏移
文章编号 0001-5733(2005)02-0415-08 中图分类号 P631 收稿日期 2004-04-13,2004-12-08收修定稿
基金项目 国家自然科学基金项目(40174017,40235055)资助.
作者简介 王红落,男,1975年生,1998年毕业于大庆石油学院,现为中国科学院地质与地球物理研究所博士研究生,主要从事勘探地震数
字处理方面的工作.E -mail:whl@https://www.wendangku.net/doc/4d4610978.html,
Receiver function forward modeling and migration based on wave -equation
finite difference method
W ANG Hong -Luo 1
,CHANG Xu 1
,C HEN Chuan -Ren
2
1Institu te o f Geolog y &Geo ph ysics ,Ch inese Aca demy o f Sc ienc es ,Bei jin g 100029,Chin a 2Th e Colleg e o f Geophysics an d Petroleum Resources ,Yan gt ze U ni ve rsit y ,Jing zhou 434023,Ch in a
Abstract The wave -equation finite difference algorithm is applied here as a solution for forward modeling and
migration of receiver function.With the help of exploration seismology concepts,the equivalent velocity and traveltime equation are established,therefore the for ward calculation and migration of receiver function are imple mented by wave -equa tion finite difference method.The numerical calculation indica tes that imaging of point diffrac tor and dome can be accomplished with high accurac y employing poststack wave -equation finite difference migration.Then the applicabilities for receiver func tion imaging be tween poststack finite difference migration and Kirchhoff migration are discussed acc ording to numerical results.Finally,the effects on velocity error for i m aging are analyzed.
Keywords Finite difference modeling,Receiver function array,Common c onverted point stack,Finite
difference migration
1 引 言
地震台站记录是震源时间函数、地下介质结构
以及仪器特性耦合作用的结果.接收函数是消除震
源效应以及仪器特性后,台站下方介质结构的径向响应.早期的接收函数方法是用单个固定台站远震体波资料来反演台站下方地壳的一维速度结
第48卷第2期2005年3月
地 球 物 理 学 报
C HINESE JOURNAL OF GEOPHYSICS
Vol.48,No.2Mar.,2005
构[1~3],后来逐渐用于调查和证明岩石圈深部间断面的存在[4~6].将不同点的一维速度结构在横向上连成剖面,就可以追踪解释地下速度界面分布[7].随着全球范围内数字地震台网的进一步加密,就有可能提取近于勘探地震2-D甚至3-D观测系统的接收函数阵列.尤其是最近流动地震台网的大范围应用,使得台站按需要布置成任意的观测系统成为可能.勘探地震在地面高密度数据采集的前提下,利用数字信号分析和波传播理论,得到2-D或3-D的地下精细结构[8].借鉴成熟的勘探地震学的方法,直接从接收函数阵列对地下界面成像,是最近地震学研究的新方向.
对于高密度空间采样接收函数阵列数据,将勘探地震中更精确的波动方程偏移成像方法推广到接收函数阵列,就具有明显的现实意义.近年来,在研究接收函数阵列PS震相运动学规律的基础上,前人已经做了很多开创性的工作.发展了对接收函数抽取共反射面元阵列并进行叠加的方法,大大提高了接收函数的信噪比[9~11].按照偏移成像的概念,要对转换界面成像,就必须把观测数据中的转换震相按波传播的相反方向聚焦到转换点上.Jones et al.[12]提出了把叠加后的接收函数进行射线路径逆向反传的技术,这种沿射线逆传的成像是一种高频近似方法,而非波动方程偏移.Ryberg et al.[13]提出了接收函数转换波与地震勘探反射波的等效时距关系,使反射地震数字处理方法和流程在接收函数上的应用成为可能.本文主要研究波动方程有限差分偏移算法在接收函数阵列上的应用.
本文考虑地震波动力学特征,引入相对传播时间的概念,采用波动方程有限差分算法正演接收函数,这样得到的接收函数保存了地震波的动力学信息,更接近于实际的接收函数波形记录.根据反射地震共中心点叠加原理,引入等效于反射地震的时距关系,在此基础上,本文对接收函数阵列实现了波动方程有限差分偏移,该方法与Kirchhoff积分偏移结果比较,成像精度明显改进.数值计算结果证明了本研究的可行性.
2基于波动方程有限差分的接收函数正演
一个地震台站接收到的三分量远震P波地震记录可表示为[3]
D V(t)=I(t)*S(t)*
E V(t),
D R(t)=I(t)*S(t)*
E R(t),
D T(t)=I(t)*S(t)*
E T(t),(1)其中D V(t),D R(t)和D T(t)分别表示地震记录中垂直分量、径向分量及切向分量,S(t)表示入射平面P波震源时间函数,I(t)表示仪器特征响应, E V(t)、E R(t)及E T(t)分别表示介质结构脉冲响应的垂直、径向及切向分量.由径向分量D R(t)、切向分量D T(t)分别与垂直分量D V(t)反褶积,获得的地球介质的径向、切向脉冲响应E R(t)、E T(t)即为通常所说的接收函数.因此,要正演接收函数,首先要正演三分量地震记录.
本文着重研究二维情况下两分量地震记录接收函数的计算.给定二维弹性介质模型,用二维弹性波动方程模拟P-SV波场,并从中提取模型表面各点波场随时间的变化,构成水平、垂直两分量地震记录.二维P-SV波在非均匀各向同性介质中的速度-应力弹性动力学方程为
Q v#x=
9R xx
9x+
9R xz
9z,
Q v#z=
9R xz
9x+
9R zz
9z,
R#xx=(K+2L)
9v x
9x+K
9v z
9z,
R#zz=K
9v x
9x+(K+2L)
9v z
9z,
R#xz=L
9v x
9z+L
9v z
9x,(2)其中(v x,v z)为质点两分量速度矢量,Q(x,z)为密度,(R xx,R xz,R zz)为应力张量,?v和?R分别表示速度及应力对时间的导数,K(x,z)和L(x,z)为拉梅系数[14].接收函数是从远震记录中得到的,在正演接收函数的计算中,考虑到远震的传播距离,可以假设远震波前进入到模型空间时是平面波,并将这一时刻作为震源的相对初始时刻,在模型空间范围内计算波场随时间的变化.
为了保证接收函数一致性并提高接收函数信噪比,最终得到波场的径向分量,需要对弹性波正演获得的两分量地震记录进行坐标旋转[15].如图1所示,把从XZ坐标系观测的水平、垂直两分量记录投影到LQ坐标系,其中L方向为直达P波的射线路径方向,Q方向与L方向正交,L与Q的分量分别对应(1)式中的D V(t)和D R(t).对各向同性介质,L方向分量主要是P波能量,Q方向分量是SV 波能量,坐标旋转其实质是一个波场分离的过程.
416地球物理学报(Chinese J.Geophys.)48卷
然后对径向分量D R (t )用震源时间函数进行反褶积,消除震源影响,根据(1)式的定义,就得到了接收函数E R (t ).但是在正演模拟中,由于震源时间函数统一用的是同样能量的Ricker 子波,所以不需要做反褶积及振幅均衡,可以认为坐标变换后的径向分量D R (t )就是接收函数
.
图1 远震平面P 波入射在绕射点产生
的绕射震相及(L ,Q )坐标系
Fig.1 Diffraction phases from a point di ffractor
and (L ,Q )coordinate
3 PS 转换波等效偏移距和等效叠加
速度
比较接收函数和地震道的异同,才有可能发现沟通二者的桥梁,建立针对接收函数等价于勘探地震道的概念和处理方法
[13]
(见图2)
.图2 转换震相的等效震源及偏移距
Fig.2 Equivalent source and offset of converted PS phase
图2所示地下一个平的速度界面,深度为Z ,界面上地层纵波速度为V P ,横波速度为V S .图2a 是勘探地震反射P 波射线路径,图2b 是远震P 波在转换界面处PS 转换波的产生和传播路径.在理解PS 转换震相和反射P 波之间关系的基础上,就可以建立等价于反射P 波的PS 震相等效震源位置和叠加速度的表达式.
如图2a,勘探地震记录中反射能量P 波最强,最成熟的方法是把一次反射P 波作为有效波,波传播的路径是从震源下行至反射点再上行回到地面检波器,是一个P 波的双程走时.震源到接收点的距离就是偏移距X .当反射界面为水平界面时,根据Snell 定律可以建立反射P 波走时T 和偏移距X 的关系
T =
2Z V P
2
+X 2
V 2P
.(3)
在X 远小于深度Z 时,也就是偏移距比较小的情况下,式(3)可以用Taylor 级数展开,得到包括偏移距高阶无穷小量的式(4),式(4)右边的第2项是反射地震多次覆盖共深度点CDP(Common Depth Point)叠加技术中正常时差NMO(Normal Move Out)
的近似值.CDP 叠加技术是反射地震最主要的处理技术之一,可以有效提高信噪比,压制多次波.
T =2Z V P +X 2
4ZV P
+O (X 4
).
(4)
对于远震台站记录,如图2b,最先看到的是初
至P 波震相,但往往无法确定其准确的震源位置和发震时间,而且初至P 波一般不反映地下界面信息,代表地下界面信息的是P 波通过界面时产生的PS 转换波.与反射P 波对比,PS 波自转换点以横波速度传播到地面台站,所以PS 震相是转换横波的单程走时.引入等效偏移距和等效叠加速度
[13]
,
可以对PS 转换波建立与反射地震时距关系同样形式的运动学方程,并根据C DP 叠加原理实现共转换点CCP(Common C onverted Point)叠加.
图2b 表示了PS 转换波和初至P 波时差T PS 之间简单的几何关系,用(5)式表示.
T PS =Z (
V -2
S -p 2
-V -2
P -p 2
),(5)
其中p 为射线参数.如果p 很小,则(5)式可以进
一步表示为
T PS =Z V -1S -V
-1P
+Z (V P -V S )p 2
2
+O (p 4).
(6)
可以认为PS 转换波和P 波到时差T PS ,就是转换波的运动学特征的充分描述,这是因为如果知道介质的V P 和V S ,就能对转换点准确成像.对比(6)式与(4)式,T PS 可以当作PS 震相的等效走时,引入等效速度V c 和等效偏移距X c ,得到
V c =2
V P V S
V P -V S
,
(7)
417
2期王红落等:基于波动方程有限差分算法的接收函数正演与偏移
X c =V P P V S V S 2pZ .(8)
等效速度和等效偏移距能够准确描述PS 震相的运动学特征,并可以按照反射地震处理技术,抽取共转换点接收函数排列,消除由射线参数不同和台站之间位置差异导致的NMO 时差,然后进行共转换点叠加.
远震接收函数一般震中距在23b 到80b 之间,如果用标准的地球模型如I ASP91(International Association of Seismology and Physics)
[16]
,则等价于P
波在Moho 界面以37b 到23b 入射,这个范围在Moho 深度处对(6)式舍去高阶小量是满足要求的.
4 波动方程有限差分偏移
偏移概念在勘探地震的初期已经被应用于实际资料,但真正基于波动方程数值方法的偏移是20世纪70年代初发展起来的
[17]
,现在波动方程偏移
成像理论和技术已经相当成熟完善[18]
.与反射地震
偏移原理相似,对接收函数偏移是对转换界面成像.其作用同样是消除转换波传播的效应,把在地表面观测到的地震记录沿转换波逆向延拓,直到转换波传播时间t =0为止,即满足成像条件.
二维声波标量方程为
92
U 9x 2+92
U 9z 2=1V 292
U
9t
2,(9)
式中U 表示波场,V 表示纵波速度,t 表示波动传
播时间,x ,z 分别表示二维空间变量.
本文在Claerbout 移动坐标变换偏移方法[19]
的
基础上对波场延拓距离进行了补偿,获得了陡倾角
界面的高精度成像[20]
.Claerbout 提出的移动坐标变
换为
x c =x ,
z c =z ,t c =
z
V
+t .(10)
式(10)代入至式(9),并考虑陡倾角偏移问题,对波场延拓距离进行补偿之后的偏移方程为
92
U c 2
+22V 92
U c 9z c 9t c =0.
(11)
反射地震叠后偏移是对(11)式求隐式有限差分解.引入等效偏移距和等效速度概念后,接收函数偏移与反射地震同样,但传播时间不是反射走时,
而是转换波与P 波时差.
5 等效叠后偏移速度
偏移成像的关键取决于速度模型的准确性.由于能够确定的只是P 波和PS 转换波的到时差,而不是PS 震相的绝对走时,这就要由反射地震偏移
的原理来确定等效偏移速度.图3是均匀介质中存在一个绕射点的情况下,在地面观测的反射波与转换波的射线路径图.左侧是双程P 波,右侧是远震入射平面P 波在经过绕射点同时产生PS 绕射转换波.
图3 均匀介质中单个绕射点零偏移距走时Fig.3 Zero offset travel time of point diffractor
i n homogeneous media
由图3的几何关系,反射地震零偏移距绕射波走时方程可表示为(12)式.对于零偏移距的接收函数,平面P 波以射线参数p =0入射,绕射S 波和直达P 波的时差如(13)式,
X 2
+Z 2
P V S 是绕射S 波
到接收点的单程走时,Z P V P 是直达P 波走时.
T P =
2V P
X 2
+Z 2
.
(12)T PS =1
V
S
X 2+Z 2-Z V P .
(13)
对式(13)加上直达P 波的时移量得到式(14),这显然是一个标准的绕射波双曲线方程,而T S 即是绕射转换波从绕射点到接收点的绝对走时.与式
(12)对比,可以得到一个等效的叠后偏移速度如(15)式.根据偏移成像原理,以V mig 作为偏移速度,做叠后时间偏移,就能使绕射转换波完全收敛到绕射点.按照(16)式进行时深转换,得到深度剖面.
T S =T PS +Z V P =1
V S X 2
+Z 2
.(14)V mig =2V S .(15)D =
V P V S
V P -V S T PS
.
(16)
418
地球物理学报(Chinese J.Geophys.)48卷
6 数值计算
首先用Virieux J.介绍的交错网格有限差分方法正演P -SV 波场的合成地震图
[14]
,边界条件用
Higdon 的吸收边界[21]
.速度模型分别模拟在相当于Moho 界面深度处存在一个绕射点和Moho 界面上的穹隆构造(图4),其目的是检验等效偏移距和等效偏移速度对接收函数成像的有效性.根据P
波速度
图4 Moho 界面的地壳速度模型:穹隆Fig.4 Crustal velocity model related to Moho
discontinuity:dome
和S 波速度的常用估计关系V P =3V S 以及Gardner
公式Q =A V B P
[22]
,由P 波速度得到S 波速度和密度.
数值模拟远震到达模型空间为平面P 波,分别
在左右两侧入射,入射角从23b 开始到37b ,每次递增1b .相当于接收到在左右两侧各15次的远震信息.震源时间函数用主频1Hz 的Ricker 子波.每个模型进行二维弹性波模拟,就可以分别得到30个XZ 两分量地震图.
利用(6)式进行CCP 叠加和偏移时,走时是PS 震相与P 波走时之差T PS ,所以要进行P 波初至拾取,把合成地震图按P 波初至排齐.然后按图1所示,需要把XZ 坐标下的波场投影到LQ 坐标系,完成直达P 波和PS 转换波波场分离,得到Q 方向上的径向分量D R (t ).从左侧30b 入射的一个记录对比来看,坐标变换后Q 方向径向分量转换波能量明显强于X 水平分量(图5).由于震源时间函数是Ricker 子波,可以不考虑震源的效应,也就是不用做反褶积来消除震源类型不同对记录的影响,也不用考虑振幅均衡问题,即D R (t )近似认为是接收函数E R (t
).
图5 平面P 波30b 入射,地下穹隆按初至排齐的水平及径向地震响应Fig.5 Seismograms of dome with first arrival as onset,generated by 30b incidence of P wave
从经过坐标投影后的D R (t )阵列中抽取共转换点接收函数(C CP)阵列,进行叠加速度分析,可以
得到与(7)式相同的叠加速度结果.再进行NMO 正常时差校正,最后叠加得到共转换点叠加接收函数阵列(图6).显然除了模型边缘覆盖次数较少,其他转换点基本都是30次覆盖,可以有效地提高信噪比.为简单起见,合成接收函数没有加随机噪声.
对比数值计算得到的共转换点叠加接收函数阵
列和速度模型,绕射点在接收函数剖面上近似一个绕射双曲线,Moho 面的突起在叠加剖面上明显比模型宽得多,这与在反射勘探地震中的现象完全吻合.从爆炸反射面假设的观点来看,叠后地震记录是地下点绕射的波场在地表叠加的结果,所以模型内各速度异常体边缘在模拟观测记录中变得模糊不清,横向分辨率降低,而偏移是波绕射的一个逆过程,使得绕射波收敛到绕射点,因此消除了由波传播效应导致的模糊化,提高横向分辨率.
419
2期王红落等:基于波动方程有限差分算法的接收函数正演与偏移
偏移速度是用模型速度换算的,即模型横波速度的2倍.图7表示了波动方程有限差分偏移方法得到的结果.与速度模型对比,成像精度很高.点绕射双曲线收敛到接近一个点,突起在横向上明显被压缩,接近模型大小,而且偏移结果的深度位置是准确的.图8是Kirchhoff积分偏移成像的结果,有限差分成像和Kirc hhoff偏移相比效果要好一些
.
图6绕射点与穹隆接收函数阵列30次覆盖共转换点叠加剖面
Fig.6Common converted point stacked seismograms of p oint diffractor and dome receiver function arrays with30
folds
图7点绕射及穹隆接收函数叠加阵列有限差分偏移成像
Fig.7Finite di fference migrations of poin t diffractor and dome stacked
arrays
图8点绕射及穹隆接收函数叠加阵列Kirchhoff偏移成像
Fig.8Kirchhoff migrations of poin t diffractor and dome stacked arrays
420地球物理学报(Chinese J.Geophys.)48卷
7 讨 论
7.1 偏移速度估计问题
对偏移成像,速度模型越准确偏移效果越好.而实际情况中速度模型往往未知.勘探地震一般从叠加速度谱,也可以辅助以声波测井资料,通过不断修正速度模型,得到最佳偏移速度.在接收函数偏移时,速度可以利用速度谱扫描分析,也可以借助区域层析成像反演的结果,但是走时层析成像方法[23]
为了能获得高精度地下结构,往往需要大规模地震数据量,其解的可靠性分析是一个困难
[24]
.速
度是影响成像效果的关键性因素,误差往往导致成
像点不能聚焦或者错误归位.图9给出了偏移速度
存在误差时的成像结果.由图所示,当偏移速度小于真实速度时,绕射双曲线未能完全收敛,造成欠偏移,成像结果是一个比绕射点浅的向上突起的Frown 波.当偏移速度大于真实速度时,绕射能量也不能在成像点很好地聚焦,而是在真实成像点的
下方形成一个下凹的Smile 波,造成过偏移[25]
.数值计算结果表明,在存在速度误差的情况下,成像结果和反射地震所得出的结论是一致的.另一方面,即使存在一定的速度误差,但偏移也能在一定程度上使绕射波收敛,提高横向分辨率,成像的结果也要远好于叠加结果
.
图9 速度误差对接收函数阵列有限差分偏移成像的影响
Fig.9 The effects of veloci ty error on finite di fference migrati on of receiver function arrays
7.2 震源及其与台站分布问题
在实际观测中,台站分布密度往往有限,震源的分布也极不均匀.这将导致接收函数的叠加次数很低,一般不利于叠前偏移.另外接收函数阵列成像,要尽可能多地收集不同方位角、不同射线参数的远震记录,以期得到不同走向、倾向及倾角的间断面的转换波信息.另外,天然地震的震源机制和类型十分复杂,在资料处理中不可能完全消除震源的效应.噪声也是不可避免的,数据本身的质量会导致速度分析的误差,并进而影响成像结果.当震源能量强弱差异比较大时,有必要进行振幅均衡.
在实际中,往往面对的是更复杂的模型,需要探测Moho 面以下更深的层位(如410km 和660km 间断面).这些都给我们用接收函数阵列偏移成像提出了更多的问题.
8 结 论
随着固定台站的加密以及流动台站的大规模使用,获取接近于勘探地震观测系统、在地表规则密集的观测数据已经成为可能,这样就可以利用接收函数阵列直接对地下速度界面成像.从勘探地震成像的原理出发,对比PS 转换震相与反射地震P 波在运动学及动力学上的特征,导出等效于勘探地震的偏移距、正常时差校正速度及叠后偏移速度的概念,从而可以实现完全利用勘探地震的数字处理技术对接收函数阵列实现偏移成像.
本文用波动方程有限差分解进行数值模拟,得到两个简单的地壳模型合成两分量记录.本文成功地利用了波动方程有限差分算法合成接收函数阵列并对其进行了偏移成像.这一方法对接收函数在空间域偏移成像的实际应用具有重要意义.
421
2期王红落等:基于波动方程有限差分算法的接收函数正演与偏移
参考文献(References)
[1]Langs ton C A.Corvallis Ore gon.crus tal and upper mantle receiver
s tructure from teles eis mic P and S waves.Bull.Seism.Soc.
Am.,1977,67:713~724
[2]Burdick L J,Langs ton C A.Modeling crusta-l s tructure through the
use of converted phases in teles eis mic body-waveforms.Bull.
Seism.Soc.Am.,1977,67:667~691
[3]Langs ton C A.Structure under Mount Rainier,Was hington,
inferred from teleseis mic body waves.J.Geophys.Res.,1979,
84(B9):4749~4762
[4]Bos tock M G.M antle s tratigraphy and evoluti on of sl ave province.
J.Geophy.Res.,1998,103(B9):21183~212000
[5]Dueker K G,Sheehan A F.Mantle discontinuity structure from
midpoi nt stacks of converted P to S waves across Yellowstone hotspot
track.J.Ge ophys.Re s.,1997,102(B4):8313~8327
[6]Kosarev G,Kind R,Sobolev S V,et al.Seis mic evidence for a
de tached Indian li thospheric mantle beneath Tibet.Science,1999,
283(5406):1306~1309
[7]吴庆举,曾融生.用宽频带远震接收函数研究青藏高原的
地壳结构.地球物理学报,1998,41(5):669~679
Wu Q J,Zeng R S.The crus tal structure of Qingha-i Xizang Plateau
inferred fro m broadband teleseis mic waveform.Chinese J.Ge ophys.
(in Chi nese),1998,41(5):669~679
[8]渥#伊尔马滋.地震数据处理.黄绪德,袁明德译.北京:石
油工业出版社,1994
Yil maz O.Seis mic Data Processi ng(in Chinese).Translated by
Huang X D,Yuan M D.Beijing:Petroleum Industry Press,1994 [9]Dueker K G,Sheehan A F.M antle di scontinui ty s tructure beneath
the Colorado rocky mountains and high plains.J.Geo phy.Res.,
1998,103(B4):7153~7169
[10]Gurrola H,Minis ter J B,Owens T.The use of velocity s pectrum for
s tacking recei ver functi on and imaging upper mantle disc onti nui ties.
Geo phys.J.Int.,1994,117(2):427~440
[11]Gurrola H,Minis ter J B.Thickness estimates of the upper mantle
transition zone from bootstrapped velocity s pectrum s tacks of receiver
functi ons.Geo phys.J.Int.,1998,133(1):31~43
[12]Jones C H,Phinney R A.Sei smic structure of the lithosphere from
teleseismic converted arrivals at s mall arrays in the southern Sierra
Nevada and vicinity,Cali fornia.J.Geo phys.Res.,1998,103
(B5):10065~10090
[13]Ryberg T,Weber M.Receiver function arrays:a reflecti on seis mic
approach.Geophys.J.Int.,2000,141(1):1~11
[14]Virieux J.P-SV wave propagati on in heterogeneous media:
Vel ocity-stress finite difference method.Geo physics,1986,51(4):
889~901
[15]Yuan X,Ni J,Kind R,et al.Lithosperic and upper mantle
structure of southern Tibet from a seis mological pas sive s ource
e xperiment.J.Geo phys.Re s.,1997,102(B12):27491~27500
[16]Kennett B L N,Engdahl E R.Traveltimes for global earthquake
l ocation and phase identi fication.Ge ophys.J.Int.,1991,105
(2):429~465
[17]Claerbout J F,D oherty S M.Downward continuati on of moveou-t
corrected seis mograms.Geophysic s,1972,37(5):741~768 [18]Gray S H,Etgen J,Dellinger J,et al.Seis mic migration problems
and s ol utions.Ge ophysics,1985,66(5):1622~1640
[19]Claerbout J F.Fundamentals of Geophysical D ata Processing.Ne w
York:Mc Graw-Hill Book Inc..1976
[20]Chang Xu,Ashida Yuz uru,Sas sa Koichi.Errors in downward
conti nuation and its correcti on in migration using a moving
coordinate transform.Geo physical Exploration(in Japanes e&
Englis h),1990,43(5):319~329
[21]Hi gdon R L.Absorbing boundary conditions for elastic waves.
Ge ophysics,1992,56(2):231~241
[22]Gardner G H F,Gardner L W,Gregory A R,et al.Formation
velocity and density-The diagnos tic basics for strati graphic traps.
Ge ophysics,1974,39(6):770~780
[23]Li u Y K,Chang X,Li u F T,et al.Three-di mensional velocity
i mages beneath the Kang-Di an Tethyan tec tonic zone of China.
Canadian Journal of Earth Scie nce s,2002,39(10):1517~1525 [24]刘伊克,常旭.地震层析成像反演中解的定量评价及应
用.地球物理学报,2000,43(2):251~256
Li u Y K,Chang X.Quanti tative assess ment of inversion s oluti on of
sei smic tomographys and its application.Chine se J.Geo phys.(i n
Chi nese),2000,43(2):251~256
[25]Zhu J,Li nes L,Gray S.Smiles and frowns in mi grati on P velocity
analysi s.Geo physics,1998,63(4):1200~1209
422地球物理学报(Chinese J.Geophys.)48卷
(完整版)大连理工大学高等数值分析抛物型方程有限差分法
抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??，T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： （1.2） ()()x x u ?=0,, ∞<<∞-x 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
现在考虑边值问题（1.1），（1.3）的差分逼近 取 N l h = 为空间步长，M T = τ为时间步长，其中N ，M 是 自然数， jh x x j ==, ()N j ,,1,0Λ=； τ k y y k ==, ()M k ,,1,0Λ= 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表 示网格节点； h G 表示网格内点（位于开矩形G 中的网格节点）的集合； h G 表示位于闭矩形G 中的网格节点的集合； h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解，N j ≤≤0，M k ≤≤0。 注意到在节点()k i t x ,处的微商和差商之间的下列关系 （(,)k j k j u u x t t t ????≡ ? ????）： ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式
常微分方程和偏微分方程的数值解法教学大纲
上海交通大学致远学院 《常微分方程和偏微分方程的数值解法》教学大纲 一、课程基本信息 课程名称（中文）：常微分方程和偏微分方程的数值解法 课程名称（英文）：Numerical Methods for Ordinary and Partial Differential Equations 课程代码：MA300 学分 / 学时：4学分 / 68学时 适用专业：致远学院与数学系相关专业 先修课程：偏微分方程，数值分析 后续课程：相关课程 开课单位：理学院数学系计算与运筹教研室 Office hours: 每周二19：00—21：00，地点：数学楼1204 二、课程性质和任务 本课程是致远学院和数学系应用数学和计算数学方向的一门重要专业基础课程，其主要任务是通过数学建模、算法设计、理论分析和上机实算“四位一体”的教学方法，使学生掌握常微分方程与偏微分方程数值解的基本方法、基本原理和基本理论，进一步提升同学们利用计算机解决实际问题的能力。在常微分方程部分，将着重介绍常微分方程初值问题的单步法，含各类Euler方法和Runge-Kutta方法，以及线性多步法。将简介常微分方程组和高阶常微分方程的数值方法。在偏微分方程部分，将系统介绍求解椭圆、双曲、抛物型方程的差分方法的构造方法和理论分析技巧，对于椭圆型方程的边值问题将介绍相应变分原理与有限元方法。将在课堂上实时演示讲授的核心算法的计算效果，以强调其直观效果与应用性。本课程重视实践环节建设，学生要做一定数量的大作业。 三、教学内容和基本要求 第一部分：常微分方程数值解法 1 引论 1.1回顾：一阶常微分方程初值问题及解的存在唯一性定理
抛物形扩散方程的有限差分法及数值实例
偏微分方程数值解 所在学院：数学与统计学院 课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：向聘
抛物形扩散方程的有限差分法及数值实例 1.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程： 22(),0u u a f x t T t x ??=+<≤?? (1.1.1) 其中a 是常数，()f x 是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, ∞<<∞-x (1.1.2) 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, 0x l << (1.1.3) 及边值条件 ()()0,,0==t l u t u ， T t ≤≤0 (1.1.4) 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 下面考虑如下热传导方程 22()(0.)(,)0(,0)()u u a f x t x u t u L t u x x ????=+????? ==??=??? (1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。 取N l h = 为空间步长，M T =τ为时间步长，其中N ，M 是自然数，用两族
平行直线jh x x j ==, ()N j ,,1,0 =和k t t k τ ==, ()M k ,,1,0 =将矩形域 G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 (),j k x t 表示网格节点；h G 表示 网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。 现在对方程进行差分近似： (一) 向前差分格式 =-+τ k j k j u u 111 2 2(())k k k j j j j j j u u u a f f f x h +--++= (1.2.2) ()j j j x u ??==0， k u 0=k N u =0 (1.2.3) 计算后得： 111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4) 其中，2 a r h τ = ，1,,1,0-=N j ，1,,1,0-=M k 。 显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下： 1000 121011000 232121000 3432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----?=+-++?=+-++??=+-++? ???=+-++? (1.2.5) 若记 () T k N k k k u u u 1 21,,,-= u ，()()()()T N x x x 121,,,-=???? ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式 10 ,0,1,,1 k k k M φ +?=+=-?=? u Au f u (1.2.6) 其中
抛物型方程的计算方法
分类号：O241.82 本科生毕业论文（设计） 题目：一类抛物型方程的计算方法 作者单位数学与信息科学学院 作者姓名 专业班级2011级数学与应用数学创新2班 指导教师 论文完成时间二〇一五年四月
一类抛物型方程的数值计算方法 (数学与信息科学学院数学与应用数学专业2011级创新2班) 指导教师 摘要: 抛物型方程数值求解常用方法有差分方法、有限元方法等。差分方法是一种对方程直接进行离散化后得到的差分计算格式，有限元方法是基于抛物型方程的变分形式给出的数值计算格式.本文首先给出抛物型方程的差分计算方法，并分析了相应差分格式的收敛性、稳定性等基本理论问题.然后，给出抛物型方程的有限元计算方法及理论分析. 关键词：差分方法，有限元方法，收敛性，稳定性 Numerical computation methods for a parabolic equation Yan qian (Class 2, Grade 2011, College of Mathematics and Information Science) Advisor: Nie hua Abstract: The common methods to solve parabolic equations include differential method, finite element method etc. The main idea of differential method is to construct differential schemes by discretizing differential equations directly. Finite element scheme is based on the variational method of parabolic equations. In this article, we give some differential schemes for a parabolic equation and analyze their convergence and stability. Moreover, the finite element method and the corresponding theoretical analysis for parabolic equation are established. Key words: differential method, finite element method, convergence, stability
常微分方程边值问题的数值解法
第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)；
(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<
时间序列分析讲义 第01章 差分方程
第一章 差分方程 差分方程是连续时间情形下微分方程的特例。差分方程及其求解是时间序列方法的基础，也是分析时间序列动态属性的基本方法。经济时间序列或者金融时间序列方法主要处理具有随机项的差分方程的求解问题，因此，确定性差分方程理论是我们首先需要了解的重要内容。 §1.1 一阶差分方程 假设利用变量t y 表示随着时间变量t 变化的某种事件的属性或者结构，则t y 便是在时间t 可以观测到的数据。假设t y 受到前期取值1-t y 和其他外生变量t w 的影响，并满足下述方程： t t t w y y ++=-110φφ (1.1) 在上述方程当中，由于t y 仅线性地依赖前一个时间间隔自身的取值1-t y ，因此称具有这种结构的方程为一阶线性差分方程。如果变量t w 是确定性变量，则此方程是确定性差分方程；如果变量t w 是随机变量，则此方程是随机差分方程。在下面的分析中，我们假设t w 是确定性变量。 例1.1 货币需求函数 假设实际货币余额、实际收入、银行储蓄利率和商业票据利率的对数变量分别表示为t m 、t I 、bt r 和ct r ，则可以估计出美国货币需求函数为： ct bt t t t r r I m m 019.0045.019.072.027.01--++=- 上述方程便是关于t m 的一阶线性差分方程。可以通过此方程的求解和结构分析，判断其他外生变量变化对货币需求的动态影响。 1.1.1 差分方程求解：递归替代法 差分方程求解就是将方程变量表示为外生变量及其初值的函数形式，可以通过以前的数据计算出方程变量的当前值。 由于方程结构对于每一个时间点都是成立的，因此可以将(1.1)表示为多个方程： 0=t ：01100w y y ++=-φφ 1=t ：10101w y y ++=φφ t t =：t t t w y y ++=-110φφ 依次进行叠代可以得到： 1011211010110101)()1()(w w y w w y y ++++=++++=--φφφφφφφφ 0111122113121102)1(w w w y y φφφφφφφ++++++=- i t i i t t i i t w y y ∑∑=-=++=0 111 1 0φφφφ (1.2) 上述表达式(1.2)便是差分方程(1.1)的解，可以通过代入方程进行验证。上述通过叠代将 t y 表示为前期变量和初始值的形式，从中可以看出t y 对这些变量取值的依赖性和动态变化 过程。 1.1. 2. 差分方程的动态分析：动态乘子(dynamic multiplier) 在差分方程的解当中，可以分析外生变量，例如0w 的变化对t 阶段以后的t y 的影响。假设初始值1-y 和t w w ,,1 不受到影响，则有：
差分方程方法
第四章 差分方程方法 在实际中，许多问题所研究的变量都是离散的形式，所建立的数学模型也是离散的，譬如，像政治、经济和社会等领域中的实际问题。有些时候，即使所建立的数学模型是连续形式，例如像常见的微分方程模型、积分方程模型等等，但是，往往都需要用计算机求数值解。这就需要将连续变量在一定条件下进行离散化，从而将连续型模型转化为离散型模型，因此，最后都归结为求解离散形式的差分方程解的问题。关于差分方程理论和求解方法在数学建模和解决实际问题的过程中起着重要作用。 下面就不同类型的差分方程进行讨论。所谓的差分方程是指：对于一个数列{}n x ，把数列中的前1+n 项()n i x i ,2,1,0=关联起来所得到的方程。 4．1常系数线性差分方程 4.1.1 常系数线性齐次差分方程 常系数线性齐次差分方程的一般形式为 02211=+?+++---k n k n n n x a x a x a x (4.1) 其中k 为差分方程的阶数，()k i a i ,,2,1 =为差分方程的系数，且()n k a k ≤≠0。对应的代数方程 02 211=++++--k k k k a a a λ λλ （4.2） 称为差分方程的（4.1）的特征方程，其特征方程的根称为特征根。 常系数线性齐次差分方程的解主要是由相应的特征根的不同情况有不同的形式。下面分别就特征根为单根、重根和复根的情况给出差分方程解的形式。 1. 特征根为单根 设差分方程（4.1）有k 个单特征根 k λλλλ,,,,321 ，则差分方程（4.1）的通解为 n k k n n n c c c x λλλ+++= 2211， 其中k c c c ,,,21 为任意常数，且当给定初始条件 () 0 i i x λ= ()k i ,,2,1 = (4.3) 时，可以唯一确定一个特解。 2. 特征根为重根 设差分方程（4.1）有l 个相异的特征根()k l l ≤≤1,,,,321λλλλ 重数分别为 l m m m ,,,21 且k m l i i =∑=1 则差分方程（4.1）的通解为
差分方程方法
第三章 差分方程方法 3.1 差分方程的平衡点及其稳定性 设有未知序列{}n x ，称 0),,,;(1=++k n n n x x x n F (3.1) 为k 阶差分方程。若有)(n x x n =，满足 0))(,),1(),(;(=++k n x n x n x n F 则称)(n x x n =是差分方程(3.1)的解，包含k 个任意常数的解称为(3.1)的通解， 110,,,-k x x x 为已知时，称其为(3.1)的初始条件，通解中的任意常数都由初始条件确定后 的解称为(3.1)的特解。 形如 )()()(11n f x n a x n a x n k k n k n =+++-++ (3.2) 的差分方程，称为k 阶线性差分方程。)(n a i 为已知系数，且0)(≠n a k 。 若差分方程(3.2)中的0)(=n f ，则称差分方程(3.2)为k 阶齐次线性差分方程，否则称为k 阶非齐次线性差分方程。 若有常数α是差分方程(3.1)的解，即0),,,;(=ααα n F ，则称α是差分方程(3.1)的平衡点，又对差分方程(3.1)的任意由初始条件确定的解)(n x x n =，都有 )(∞→→n x n α，则称这个平衡点α是稳定的。 若110,,,-k x x x 已知，则形如),,,;(11-+++=k n n n k n x x x n g x 的差分方程的解可以在计算机上实现。下面给出理论上需要的一些特殊差分方程的解。 一阶常系数线性差分方程 b x x n n =++α1， (3.3) （其中b ,α为常数，且0,1-≠α）的通解为 )1()(++-=a b C x n n α (3.4) 易知)1(+αb 是方程(3.3)的平衡点，由(3.4)式知，当且仅当1<α时，)1(+αb 是稳定的平衡点。
差分方程模型的理论和方法
第九章 差分方程模型的理论和方法 引言 1、差分方程： 差分方程反映的是关于离散变量的取值与变化规律。通过建立一个或几个离散变量取值所满足的平衡关系，从而建立差分方程。 差分方程就是针对要解决的目标，引入系统或过程中的离散变量，根据实际背景的规律、性质、平衡关系，建立离散变量所满足的平衡关系等式，从而建立差分方程。通过求出和分析方程的解，或者分析得到方程解的 特别性质（平衡性、稳定性、渐近性、振动性、周期性等），从而把握这个离散变量的变化过程的规律，进一步再结合其他分析，得到原问题的解。 2、应用：差分方程模型有着广泛的应用。实际上，连续变量可以用离散变量来近似和逼近，从而微分方程模型就可以近似于某个差分方程模型。差分方程模型有着非常广泛的实际背景。在经济金融保险领域、生物种群的数量结构规律分析、疾病和病虫害的控制与防治、遗传规律的研究等许许多多的方面都有着非常重要的作用。可以这样讲，只要牵涉到关于变量的规律、性质，就可以适当地用差分方程模型来表现与分析求解。 3、差分方程建模： 在实际建立差分方程模型时，往往要将变化过程进行划分，划分成若干时段，根据要解决问题的目标，对每个时段引入相应的变量或向量，然后通过适当假设，根据事物系统的实际变化规律和数量相互关系，建立每两个相邻时段或几个相邻时段或者相隔某几个时段的量之间的变化规律和运算关系（即用相应设定的变量进行四则运算或基本初等函数运算或取最运算等）等式（可以多个并且应当充分全面反映所有可能的关系），从而 建立起差分方程。或者对事物系统进行划分，划分成若干子系统，在每个子系统中引入恰当的变量或向量，然后分析建立起子过程间的这种量的关系等式，从而建立起差分方程。在这里，过程时段或子系统的划分方式是非常非常重要的，应当结合已有的信息和分析条件，从多种可选方式中挑选易于分析、针对性强的划分，同时，对划分后的时段或子过程，引入哪些变量或向量都是至关重要的，要仔细分析、选择，尽量扩大对过程或系统的数量感知范围，包括对已有的、已知的若干量进行结合运算、取最运算等处理方式，目的是建立起简洁、深刻、易于求解分析的差分方程。在后面我们所举的实际例子中，这方面的内容应当重点体会。 差分方程模型作为一种重要的数学模型，对它的应用也应当遵从一般的数学建模的理论与方法原则。同时注意与其它数学模型方法结合起来使用，因为一方面建立差分方程模型所用的数量、等式关系的建立都需要其他的数学分析方式来进行；另一方面，由差分方程获得的结果有可以进一步进行优化分析、满意度分析、分类分析、相关分析等等。 第一节 差分方程的基本知识 一、 基本概念 1、 差分算子 设数列{}n x ，定义差分算子n n n x x x -=??+1:为n x 在n 处的向前差分。 而1--=?n n n x x x 为n x 在n 处的向后差分。 以后我们都是指向前差分。 可见n x ?是n 的函数。从而可以进一步定义n x ?的差分： n n x x 2)(?=?? 称之为在n 处的二阶差分，它反映的是的增量的增量。 类似可定义在n 处的k 阶差分为：
双曲方程基于matlab的数值解法
双曲型方程基于MATLAB 的数值解法 （数学1201,陈晓云,41262022） 一：一阶双曲型微分方程的初边值问题 0,01,0 1.(,0)cos(),0 1. (0,)(1,)cos(),0 1. u u x t t x u x x x u t u t t t ππ??-=≤≤≤≤??=≤≤=-=≤≤ 精确解为 ()t x cos +π 二：数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题，将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ为空间和时 间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2：差分格式的建立 0u u t x ??-=?? 2.2.1：Lax-Friedrichs 方法 对时间、空间采用中心差分使得 2h 1 1111)(2 1u u x u u u u u t u k j k j k j k j k j k j -+-++-= +=-= ????τ τ 则由上式得到Lax-Friedrichs 格式 1 11111()202k k k k k j j j j j u u u u u h τ+-+-+-+-+=
截断误差为 ()[]k k k j h j j R u L u Lu =- 1 11111()22k k k k k k k j j j j j j j u u u u u u u h t x τ+-+-+-+-??=+-+?? 23222 3 (),(0,0)26k k j j u u h O h j m k n t x ττ??= -=+≤≤≤≤?? 所以Lax-Friedrichs 格式的截断误差的阶式2()O h τ+ 令/s h τ=：则可得差分格式为 1111 11(),(0,0)222 k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤ 0cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤ 其传播因子为： ()()()e e G h i h i s h i h i σσσστσ---=-+e e 221， 化简可得： ()()()()()h s G h is h G στσσστ σsin 11,sin cos ,2 2 2--=-= 所以当1s ≤时，()1,≤τσG ,格式稳定。 * 2.2.2：LaxWendroff 方法 用牛顿二次插值公式可以得到LaxWendroff 的差分格式，在此不详细分析，它的截断误差为() h 2 2 +O τ ，是二阶精度；当2s ≤时，()1,≤τσG ， 格式稳定。在这里主要用它与上面一阶精度的Lax-Friedrichs 方法进行简单对比。 2.3差分格式的求解
第九章-偏微分方程差分方法汇总
第9章 偏微分方程的差分方法 含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。 9.1椭圆型方程边值问题的差分方法 9.1.1 差分方程的建立 最典型的椭圆型方程是Poisson （泊松）方程 G y x y x f y u x u u ∈=??+??-≡?-),(),,()(2222 （9.1） G 是x ,y 平面上的有界区域，其边界Γ为分段光滑的闭曲线。当f (x ,y )≡0时，方程 （9.1）称为Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件 第一边值条件 ),(y x u α=Γ （9.2） 第二边值条件 ),(y x n u β=??Γ （9.3） 第三边值条件 ),()( y x ku n u γ=+??Γ （9.4） 这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数，k (x,y )≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。 用差分方法求解偏微分方程，就是要求出精确解u (x ,y )在区域G 的一些离散节点（x i ，y i ）上的近似值u i ,j ≈(x i ，y i ）。差分方法的基本思想是，对求解区域G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。 设G ={0 i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 双曲型方程的有限差分法 线性双曲型方程定解问题： （a ）一阶线性双曲型方程 ()0=??+??x u x a t u （b ）一阶常系数线性双曲型方程组 0=??+??x t u A u 其中A ，s 阶常数方程方阵，u 为未知向量函数。 （c ）二阶线性双曲型方程（波动方程） ()022=?? ? ??????-??x u x a x t u ()x a 为非负函数 （d ）二维，三维空间变量的波动方程 0222222=???? ????+??-??y u x u t u 022222222=???? ????+??+??-??z u y u x u t u §1 波动方程的差分逼近 1.1 波动方程及其特征 线性双曲型偏微方程的最简单模型是一维波动方程： （1.1） 22 222x u a t u ??=?? 其中0>a 是常数。 （1.1）可表示为：022 222=??-??x u a t u ，进一步有 0=??? ????+?????? ????-?? u x a t x a t 由于 x a t ?? ±??当a dt dx ±=时为()t x u ,的全导数 （=dt du dt dx x u t u ???+??x u a t u ??±??=）,故由此定出两个方向 （1.3） a dx dt 1±= 解常微分方程（1.3）得到两族直线 （1．4） 1C t a x =?+ 和 2C t a x =?- 称其为特征。 特征在研究波动方程的各种定解问题时，起着非常重要的作用。 比如，我们可通过特征给出（1.1）的通解。（行波法、特征线法） 将（1.4）视为),(t x 与),(21C C 之间的变量替换。由复合函数的微分法则 2 12211C u C u x C C u x C C u x u ??+ ??=?????+?????=?? x C C u C u C x C C u C u C x u ????? ? ????+????+?????? ????+????=??2 212121122 2 22122212212C u C C u C C u C u ??+???+???+??= 22 22122122C u C C u C u ??+???+??= 同理可得 a t t a t C -=??-=??1，a t C =??2 ???? ????-??=?????+?????=??21 2211C u C u a t C C u t C C u t u 毕业论文 题目抛物型方程的差分解法学院数学科学学院 专业信息与计算科学 班级计算0802 学生王丹丹 学号20080901045 指导教师王宣欣 二〇一二年五月二十五日 摘要 偏微分方程的数值解法在数值分析中占有重要的地位，很多科学技术问题的数值计算包括了偏微分方程的数值解问题【1】。近三十多年来，数值解法的理论和方法都有了很大的发展，而且在各个科学技术的领域中应用也愈来愈广泛。本文的研究主要集中在依赖于时间的问题，借助于简单的常系数扩散方程，介绍抛物型方程的差分解法。本文以基本概念和基本方法为主，同时结合算例实现算法。 第一部分介绍偏微分方程及差分解法的基本概念，引入本文的研究对象——常系 数扩散方程： 2 2 ,,0 u u a x R t t x ?? =∈>?? 第二部分介绍上述方程的几种差分格式及每种格式的相容性、收敛性与稳定性。 第三部分通过算例检验每种差分格式的可行性。 关键词：偏微分方程；抛物型；差分格式；收敛性；稳定性；算例 ABSTRACT The numerical solution of partial differential equation holds an important role in numerical analysis .Many problems of compution in the field of science and techology include the numerical solution of partial differential equation. For more than 30 years, the theory and method of the numerical computation made a great development and its applications in various fields of science and technology are more and more widely. This paper focuses on the problems based on time. I will use object-constant diffusion equation to introduces the finite difference method of parabolic equation. This paper mainly focus on the basic concept ,basic method and simple numerical example. The first part of this paper introduces partial differential equations and basic concepts of finite difference method.I will introduce the object-constant diffusion equation for the first time. 2 2 ,,0 u u a x R t t x ?? =∈>?? The second part of this paper introduces several difference schemes of the above equation and their compatibility ,convergence and stability. The third part tests the accuracy of each scheme. Key words:partial differential equation；parabolic；difference scheme；convergence；stability；application 第三节 差分方程常用解法与性质分析 1、常系数线性差分方程的解 方程)(...110n b x a x a x a n k k n k n =+++-++ （ 8） 其中k a a a ,...,,10为常数，称方程（8）为常系数线性方程。 又称方程0...110=+++-++n k k n k n x a x a x a （9） 为方程（8）对应的齐次方程。 如果（9）有形如 n n x λ=的解，带入方程中可得： 0 ...1110=++++--k k k k a a a a λλλ （10） 称方程（10）为方程（8）、（9）的特征方程。 显然，如果能求出（10）的根，则可以得到（9）的解。 基本结果如下： （1） 若（10）有k 个不同的实根，则（9）有通解： n k k n n n c c c x λλλ+++=...2211， （2） 若（10）有m 重根λ，则通解中有构成项： n m m n c n c c λ )...(121----+++ （3）若（10）有一对单复根 βαλi ±=，令：?ρλi e ±=， αβ?βαρarctan ,22=+=，则（9）的通解中有构成项： n c n c n n ?ρ?ρsin cos 21--+ （4） 若有m 重复根：βαλi ±=，φρλi e ±=，则（9）的通项中有成 项： n n c n c c n n c n c c n m m m m n m m ?ρ?ρsin )...(cos )...(1221121---++---+++++++ 综上所述，由于方程（10）恰有k 个根，从而构成方程 （9）的通解中必有k 个独立的任意常数。通解可记为：-n x 如果能得到方程（8）的一个特解：*n x ，则（8）必有通解： =n x -n x +* n x （11） （1） 的特解可通过待定系数法来确定。 例如：如果)(),()(n p n p b n b m m n =为n 的多项式，则当b 不是特征 根时，可设成形如)(n q b m n 形式的特解，其中)(n q m 为m 次多项式；如 果b 是r 重根时，可设特解：r n n b )(n q m ，将其代入（8）中确定出系 数即可。 第五章 常微分方程的差分方法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：欧拉公式、改进的欧拉公式。 三、教学重点难点 1．教学重点：改进的欧拉公式、龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 基于数值积分的求解公式：欧拉公式、改进的欧拉公式 引 言 1．主要考虑如下的一阶常微分方程初值问题的求解: 00()(,)()y x f x y y x y '=??=? 微分方程的解就是求一个函数y=y(x)，该函数满足微分方程并且符合初值条件。 2. 例如微分方程: xy'-2y=4x ；初始条件: y(1)=-3。 于是可得一阶常微分方程的初始问题 24(1)3y y x y ?'=+???=-?。 显然函数y(x)=x 2-4x 满足以上条件,因而是该初始问题的微分方程的解。 3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。因此，只能依赖于数值方法去获得微分方程的数值解。 4． 微分方程的数值解： 设微分方程问题的解y(x)的存在区间是[a,b]，初始点x 0=a ，将[a,b]进行划分得一系列节点x 0 , x 1 ,...,x n ，其中a= x 0< x 1<…< x n =b 。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点x k 的近似值y(x k )，即 y≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 如果计算y n 时，只利用y n-1，称这种方法为单步法；如果在计算y n 时不仅利用y n-1，而且还要利用y n-2, y n-3,…, y n-r ,则称这种方法为r 步方法，也称多步法。 §5.1 欧拉方法 §5.1.1 欧拉格式 方程()(,)n n n y x f x y '=中，1()()()n n n y x y x y x h +-'≈ 1()()(,())n n n n y x y x hf x y x +≈+?1(,)n n n n y y hf x y +=+ 称为解一阶常微分方程初值问题的欧拉公式，也称显示欧拉公式。 欧拉公式的几何意义非常明显，因为微分方程的解在xoy 平面上表示一族积分曲线。用欧拉公式求数值解的几何意义如图： 容易验证，该折线各个顶点的纵坐标(1,2...)n y n =就是欧拉公式算得的近似值解，所以，欧拉方法又称为折线法。 算例：P98常微分方程数值解法
大连理工大学 高等数值分析 偏微分方程数值解(双曲方程书稿)
常微分方程差分解法、入门、多解法
差分方程的解法
常微分方程的差分方法