v2第9章常微分方程数值解wb.docx
- 文档编号:13165880
- 上传时间:2023-06-11
- 格式:DOCX
- 页数:70
- 大小:180.88KB
v2第9章常微分方程数值解wb.docx
《v2第9章常微分方程数值解wb.docx》由会员分享,可在线阅读,更多相关《v2第9章常微分方程数值解wb.docx(70页珍藏版)》请在冰点文库上搜索。
v2第9章常微分方程数值解wb
第九章常微分方程数值解法
在工程和科学技术的实际问题中,常需要求解常微分方程,许多实际问题的数学模型都是常微分方程
或常微分方程的定解问题,如物体运动、电路振荡、化学反映及生物群体的变化等。
例如:
对于单摆问题
(如图9.1),选取摆长I
的函数,来描述单摆运动。
到下面二阶常微分方程
1和质量m1,我们希望得到摆角q的关于时间t
根据力矩与角加速度的关系,
经一定的简化后可得
d2qdt2
gsinq
其中g为重力加速度。
当单摆在摆动开始时刻tt0时的初始摆角qt0
qo
和初始角速度凹
dttto
q(t0)q0都确定时,
单摆运动规律q(t)惟一确定,
这就可以写成一个初值问题:
d2q
dF
gsinq
qtoqo,q(to)q。
虽然常微分方程广泛存在于众多研究和应用领域,但是常微分方程中往往只有少数较简单和典型的微
分方程(例如线性常系数常微分方程等)的解能用初等函数、特殊函数或它们的级数与积分表达,对于变系
数常微分方程的解析求解就比较困难,而一般的非线性常微分方程的求解困难就更不用说了。
大多数情况下,常微分方程只能用近似方法求解。
这种近似解法可分为两大类:
一类是近似解析法,如级数解法、逐
次逼近法等;另一类是数值解法,它给出方程在一些离散点上的近似解。
在具体求解常微分方程时,需具备某种定解条件,常微分方程和定解条件合在一起组成定解问题。
定
解条件有两种:
一种是给出积分曲线在初始点的状态,称为初始条件,相应的定解问题称为初值问题;另一类是给出积分曲线首尾两端的状态,称为边界条件,相应的定解问题称为边值问题。
本章主要介绍常微分方程初值问题的基本数值方法、理论和算法。
9.1基本概念
9.1.1常微分方程初值问题
常微分方程可分为线性、非线性、高阶方程与方程组等类,其中线性方程包含于非线性类中,高阶方程可化为一阶方程组,若方程组中的所有未知量看作一个向量,则方程组可写成向量形式的单个方程。
因
此研究一阶常微分方程的初值问题
dy
f(x,y),
dx
y(a)
yo
的数值解法具有典型性,其中方程的解为
y:
R
R。
axb
(9-1)
只有保证问题(9-1)的解存在惟一的前提下,其数值解法的研究才有意义。
由常微分方程的基本理论,我们有:
定理9.1如果(9-1)中的f(x,y)满足条件
(1)f(x,y)在区域D{(x,y)axb,y}上连续;
(2)f(x,y)在D上关于y满足Lipschitz条件,即存在常数L,使得
f(x,y)f(x,y)Lyy|
则初值问题(9-1)在区间[a,b]上存在惟一的连续解yy(x)。
在本章的讨论中,我们总假定方程满足以上两个条件,从而方程总存在惟一的连续解。
在实际问题的研究中,许多问题最后可以归结为一阶常微分方程组的初值问题Yifi(x,Y1,Y2丄,Ym)
(i1,L,m)
Yi(a)Yio
若记Y(Y1,Y2,L,Ym)T,Yo(Y10,Y20,L,ymo)T,f(f1,f2,L,fm)T,则初值问题(9-2)可写成如下向量
(9-2)
形式的单个方程
如果向量函数f(x,y)在区域D:
axb,y使得对x[a,b],Y1,Y2Rm,都有
If(x,yj
那么问题(9-3)在[a,b]是存在惟一解y二y(x)。
yf(x,y)
y(a)y。
Rm连续,且关于y满足Lipschitz条件,即存在L0,
f(x,y2)L||y1y2,
问题(9-3)与问题(9-1)形式上完全相同,故对初值问题(9-1)所建立的各种数值解法都可应用于求解问题
(9-3),只需将y换成向量y,f(x,y)换成f(x,y)即可。
对于高价常微分方程的初值问题的处理,
下:
设有m阶常微分方程初值问题
可以通过变量代换化为一阶常微分方程组初值问题,过程如
(9-3)
引入新变量y1
y,y2
y,l,ym
(m)
y
y(a)
y(m"
y1
y2
y2
y3
M
(m1)、
f(x,Y,Y,L,Y),axb
(1)(m1)(m1)
Yo,Y(a)Yo,L,Y(a)Yo
,问题(9-4)就化为一阶常微分方程组初值问题
%(a)Yo
Y2(a)y0°
M
Ym1(a)y0m2)
Ym(a)y0m1}
ym1Ym
Ymf(x,Y1,Y2,L
ym)
(9-4)
令
y(Y1,Y2,l,ym1,ym)T,Yo(Yo,y01},l,y0m
2)
m1)、T
T
)',f(y2,yo,L,Ym,f(x,儿y2,L,ym))
则该一阶微分方程组初值问题可写成(9-3)的向量形式。
对于问题(9-1),在区间a,b
9.1.2初值问题数值解法
所谓数值解法,是通过常微分方程离散化而给出解在某些节点上的近似值。
引入若干节点
axoxiX2LXnb
hnXn1Xn,(n0,1丄,N1)称为由Xn到Xni的步长,当hnh(常数)时称为定步长,否则称为变
ba
步长。
多数情况下,采用等步长,即h,Xnanh(n1,2,L,N)。
记问题(9-1)的精确解为y(x),
N
记y(Xn)的近似值为yn,记f(Xn,yn)为fn,Yn51,2,,N)称为问题(9-1)的数值解。
求初值问题数值解的方法一般采用步进法,即在计算出y,in后计算yn1,方法分为单步法和多步法。
单步法是指在计算yn1时只利用yn,而多步法在计算yn1时不仅要利用y.,还要利用前面已经计算出来的
若干个ynj,j1,2丄,11。
我们称要用到yn和ynj,j1,2,L,l1的多步方法为I步方法。
不论单步法和多步法,它们都有显式方法和隐式方法之分。
显式单步法的计算公式可以写为
yn+1=ynh(Xn,yn,h)(9-5)
此公式右端不含yn1。
隐式单步法的计算公式可以写为
yn+1二ynh(Xn,yn,yn1,h)(9-6)
在(9-6)中右端含有yn1,从而(9-6)是yn1的方程式,需要求解方程或者采用迭代法。
显式多步法的计算公式可以写为
yn+1二ynh(Xn,yn,yn-1,L,y.l1,h)
隐式多步法的计算公式可以写为
yn+1二ynh(Xn,yn1,yn,yn-1,L,ynI1,h),kI1
显式公式与隐式公式各有特点。
显式公式的优点是使用方便,计算简单,效率高,其缺点是计算精度
低,稳定性差。
隐式公式正好与它相反,它具有计算精度高,稳定性好等优点,但求解过程很复杂,一般
采用迭代法。
为了结合各自的优点,通常将显式公式与隐式公式配合使用,由显式公式提供迭代初值,再经隐式公式迭代校正。
9.2欧拉方法
欧拉公式是常微分方程初值问题数值解法中最简单的方法,它的导出过程能较清楚地说明常微分方程建立数值解法的基本思想,下面我们介绍欧拉方法。
9.2.1欧拉方法
欧拉方法的建立可以通过下面的三种方法:
、用差商近似导数
在问题(9-1)中,若用向前差商也Q代替y(Xn),则得
h
y(Xn1)y(Xn)
h
f(Xn,yn(Xn))
(n
y(Xn)用其近似值yn代替,所得结果作为
y(Xn1)的近似值,记为
0,1,,N1)
yn1,则有
yn1ynhf(Xn,yn)(n0,1,L,N1)
这样,问题(9-1)的近似解可以表示为
£
5
—
A
—
L弋
旳Xj0
9.2.2隐式欧拉方法和梯形方法
图9.2
在微分方程离散化时,如果用向后差商代替导数,即
y(Xn1)
y(Xn1)y(Xn)
则得到如下差分方程
y1ynhf(Xn,yn)(n0,1L,N1)
(9-7)yoy(xo)
按式(9-7)由初值y0经过N步迭代,可逐次算出y1,y2,yN,此方程称为差分方程。
式(9-7)称为求解一阶
常微分方程初值问题(9-1)的欧拉公式,也称显式欧拉公式。
需要说明的是,用不同的差商近似导数,将得到不同的计算公式。
二、用Taylor多项式近似
把y(Xn1)在Xn点处Taylor展开,取一次多项式近似,则得
h2
y(Xn1)y(Xn)hy(Xn)可y
y(Xn)hf(Xn,y(Xn))
h2Ny()[Xn,Xn1]
设步长h的值较小,略去余项,并以yn代替y(xj,便得
yn1ynhf(Xn,yn)
三、用数值积分法
将问题(9-1)中的微分方程在区间[xn,xn1]上两边积分,可得
Xn1
y(Xn1)y(Xn)Xf(x,y(x))dx(n0,1,,N1)(9-8)
Xn
用yn1,yn分别代替y(Xn1),y(Xn),若对右端积分采用取左端点的矩形公式,即
xn1
f(x,y(x))dxhf(Xn,yn)
Xn
同样可得到显式欧拉公式(9-7)。
类似地,对右端积分采用其它数值积分方法,又可得到不同的计算公式。
以上三种方法都是将微分方程离散化的常用方法,每一类方法又可导出不同形式的计算公式,其思想
在研究常微分方程数值解法中具有重要意义。
对于欧拉方法,从几何上看,微分方程yf(x,y)在xoy平面上确定了一个向量场:
点(x,y)处的方
向斜率为f(x,y)。
如图2.2,问题(9-1)的解yy(x)代表一条过点(x0,y0)的曲线,称为积分曲线,且此曲线上每点的切向都与向量场在这点的方向一致。
从点P0(x),y0)出发,以f(x0,y0)为斜率作一直线段,
与直线xX1交于点只(为,丫1),显然有yyhf(x°,y0),再从R出发,以彳(花,%)为斜率作直线段推进到xX2上一点P2(X2,y2),其余类推,这样得到解曲线的一条近似曲线,它就是折线PPP2L。
因
此欧拉方法又称为欧拉折线法。
上面我们给出了求解一阶常微分方程初值问题(9-1)的一种
最简单的数值公式:
欧拉公式(9-7)。
虽然它的精度比较低,实践
中很少采用,但它的导出过程能较清楚地说明构造数值解公式的基本思想,且几何意义明确,因此它在理论上仍占有一定的地位。
yiyhf(xni,yni)(n0,1,L,N1)
yoy(x。
)
(9-9)
此公式称为求解问题(9-1)所得的数值解称为隐式欧拉法。
隐式欧拉法与欧拉公式(9-7)形式上相似,但实际计算时却复杂得多。
欧拉公式(9-7)计算
不含yn1,但是隐式欧拉法计算yn1的公式中含有yn1。
在求解yn1时,yn,Xn1为已知,
yn1ynhf(Xn1,yn1)的根。
一般说来,这是一个非线性方程,当不能精确求解得到yn1
我们需要运用简单迭代法来求解。
迭代格式为
yV]1ynhf(Xn,yn)
ynk11]ynhf(Xn1,ynk]1)(k0,1,2,L)
由于f(x,y)满足Lipschitz条件,所以
yn1yn1hf(Xn1,yn1)f(Xn1,Yn1)hLYn1Yn1
由此可知,只要0hL1,迭代法就收敛到解yn1。
yn1的公式中
yn1是方程
的表达式时,
(9-8)中右端积分,即
利用数值积分方法将微分方程离散化时,若用梯形公式计算式
冷1h
xf(X,y(X))dX-[f(Xn,y(Xn))f(Xn1,y(Xn1))]
n2
并用yn,yn1代替y(Xn),y(Xn1),则得计算公式
h
yn1yn2【f(Xn,yn)f(Xn1,『n
1)]
(9-10)
这就是求解初值问题(9-1)的梯形公式。
梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为
yn0]1ynhf(Xn,yn)
f(Xn,yn)f(Xn1^^)
(k
0,1L)
由于函数f(x,y)关于y满足Lipschitz条件,所以
h
2
ynk11]yn1
hL
[k]
yn1yn1
2
yn1。
[k]
f(Xn1,yn1)f(Xn1,『n1)
hI
0——1时,迭代法就收敛到解
2
9.2.3局部截断误差与方法的精度为了刻画数值解的准确程度,引入局部截断误差与方法精度的概念。
定义9.1假设在某一步的数值解是准确的,即
用某公式推算所得yn1,我们称
其中L为Lipschitz常数。
因此,当
yny(Xn)(这个假设称为局部化假设)。
在此前提下,
Rn1
y(Xn1)yn1
为该公式(即该方法)的局部截断误差。
定义9.2如果某种方法的局部截断误差满足
Rn1y(Xn
1)yn1O(hp1)
则称该方法是p阶方法,或具有p阶精度。
显然p越大,方法的精度越高。
、欧拉法的截断误差
假设问题的解y(x)充分光滑,且前n步计算结果是精确的,即
yiy(Xi),
y(x)f(Xi,y(x))(in)
于是欧拉法的截断误差是
R1y(Xn
y(Xn
h2
i)
1)
yn1y(Xn1)ynhf(Xn,yn)y(Xn)hy(Xn)
(Xn)
O(h3)
h2
这里y(xn)称为局部截断误差主项。
显然
2
二、隐式欧拉法的截断误差计算如下:
Rn1y(Xn1)yn1丫区1)Yn
y(Xn1)y(Xn)hy(Xn+J
h23
y(Xn)+hy(Xn)+y(Xn)0(h)
2
2
Rn10(h),所以欧拉法是一阶方法。
hf(Xn+i,yn+i)
y(Xn)h(y(Xn)+hy(Xn)+^2y(Xn)0(『))
2
h^y(Xn)O(h3)
2
所以隐式欧拉法是一阶方法,其局部截断误差主项是
h2
7y(Xn)°
三、梯形法的截断误差是
Rn1
y(Xn1)yn1y(Xn1)
yn
h
2[f(Xn,yn)f(Xn1,yn1)]
h
y(Xnh)亦)门丫区)
y(Xnh)]
12y(Xn)0(h4)
12
所以梯形法是二阶方法,其局部截断误差主项是
h3
12y
(Xn)°
924改进的欧拉法
通过上面的分析,我们看到,相比于显式欧拉法和隐式欧拉法,梯形方法提高了精度,但是梯形方法,
在应用迭代公式进行实际计算时,每迭代一次,都要重新计算函数f(x,y)的值,而迭代又要反复进行若干
次,计算量很大,当函数f(x,y)比较复杂时,这个问题会变得更加突出。
为了控制计算量,我们先用欧拉
公式求得一个初步的近似值yn1,称之为预测值,预测值yn1的精度可能很差,再用梯形公式将它校正一
次得yn1,称为校正值。
这样的预测校正系统通常称为改进的欧拉公式。
即
预测:
yn1ynhf(Xn,yn)
h—
校正:
yn1yn-f(Xn,yn)f(Xn1,、n1)
为了便于编制程序上机,将上式改写成
ypynhf(Xn』n)
(9-11)
yqynhf(Xnh,yp)
yn12(ypyq)
算法实现如表9.1o
算法9.1
①
输入a,b,f(x,y),h,初值y0
②
Nb
a门
n0,xa,yyh
③
计算
Yp
yhf(x,y)
Yq
yhf(xh,Yp)
表9.1
1
-(ypyq)y,xhx
输出(X,y)
④若nN1,置n1n,转③;否则停机。
F面计算改进欧拉法的截断误差。
改进欧拉法可以改写成
h
Ya1Ya[f(Xn,yn)f(xn1,Yahf(Xn,yn))]
2
展开式,注意到YnY(Xn),有
hh
Y(Xn)-f(Xn,Y(Xn))-[f(Xn,y(Xn))
22
2
hfx(Xn,Y(Xn))hfy(Xn,y(Xn))f(Xn,y(Xn))O(h)]
h23
Y(Xn)hy(Xn)Y(Xn)O(h)
2
在(xn,y(xn))处作Taylor
Yn
3
Rn1y(Xn1)Ya1O(h)
所以改进欧拉法是二阶方法。
例9.1用欧拉法、隐式欧拉法、
梯形法和改进欧拉法解初值问题
yyx1,0x1dx
y(0)1
此问题的精确解为Y(x)e-X+xo
解取步长h0.1,对于此问题,欧拉法为
9
n10
yn+1yn
10
100
隐式欧拉法为
10
n11
yn+1yn
11
110
梯形法为
19
n1
yn+1yn
21
10510
改进欧拉法为
yn+10.905yn
0.0095n0.1
计算每种方法所得与精确解之间的误差,见表9.2。
表9.2
Xn
欧拉法
隐式欧拉法
梯形法
改进欧拉法
(Yny(Xn)|)
(
yy(Xn)|)
(
yy(Xn)|)
(yny(Xn)|)
0.0
0.0
0.0
0.0
0.0
0.1
3
4.810-
3
4.310-
-5
7.510
1.610-4
0.2
8.710-3
7.710-3
1.410-4
2.910-4
0.3
-2
1.210
-2
1.010
1.910-4
4.010-4
0.4
1.410-2
1.310-2
2.210-4
4.810-4
0.5
1.610-2
1.410-2
-4
2.510
-4
5.510
0.6
-2
1.710
-2
1.610
2.710-4
5.910-4
0.7
-2
1.8102
-2
1.7102
2.910-4
6.210-4
0.8
1.910-2
1.710-2
3.010-4
6.510-4
0.9
-2
1.910
-2
1.810
3.110-4
6.610-4
1.0
1.910-2
1.810-2
3.110-4
6.610-4
例9.2设位于坐标原点的甲舰向位于x轴上的点A(1,0)处的乙舰发射导弹,导弹始终对准乙舰。
如果乙舰
以最大的速度o(o是常数)沿平行于y轴的直线行驶,导弹的速度是5°,通过欧拉法和改进欧拉法,模
拟出导弹运行的轨道曲线,并与其精确轨道进行比较。
解设导弹的轨迹曲线为yy(x),并设经过时间t,导弹位于点P(x,y),乙舰位于点Q(1,ot)。
由于导弹始终对准乙舰。
故此时直线PQ就是导弹的运动轨迹曲线0P在点P处的切线,即有
dy°ty
dx1x
亦即
ot(1x)y'y
又根据题意,弧OP的长度为
AQ的5倍,即
y'2dx
5°t
由此得
1X2
(1x)y'y50Jy'dx
等式两边关于x求导得
(1x)y''-^1y'2
5
结合初值条件y(0)0,y'(0)0,利用常微分方程的知识,可求出此问题的精确解
5(1x)512(1
812
5
24
F面我们考虑用数值算法来求解此问题。
引入y1y,y2y,上述二阶常微分方程可以转化成一阶常微分
方程组
Y1
y2
y2,
2
y2
5(1x)
y1(0)0
y2(0)0
9.3
用欧拉法和改进欧拉法在区间0,1求解此问题,选取h0.01,方法求得的图像以及精确解的图像如图
所示:
图9.3
由图可知,欧拉法和改进欧拉法都能比较好地模拟导弹的运行轨道,并且改进欧拉法更加精确。
人物介绍
莱昂哈德欧拉(LeonhardEuler,1707-1783),瑞士数学家、自然科学家。
1707年4月15日出生于瑞士的巴塞尔,1783年9月18日于俄
国圣彼得堡去逝。
欧拉出生于牧师家庭,自幼受父亲的教育,13岁时入读巴塞尔大学,15岁大学毕业,16岁获得硕士学位。
欧拉是18世纪数学界最
杰出的人物之一,他不但为数学界做出巨大贡献,更把数学推至几乎整个物理的领域。
他是数学史上最多产的数学家,平均每年写出八百多页的论文,还写了大量的力学、分析学、几何学、变分法等课本,其中《无穷小分析引论》、《微分学原理》、《积分学原理》等都成为数学中的经典著作。
欧拉对数学的研究如此广泛,因此在许多数学的分支中也可经常见到以他的名字命名的重要常数、公式和定理。
9.3龙格-库塔(Runge-Kutta)法
在数值分析中,龙格-库塔法(Runge-Kutta)是用于求解常微分方程重要的一类隐式或显式方法,
这些技术由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明。
本节我们来学习龙格-库塔法。
9.3.1Runge-Kutta法的一般形式
设初值问题(9-1)的解yy(x)C1[a,b],由微分中值定理,我们知道,必存在[Xn’Xm],使
y(xn1)y(xn)hy()
ynhf(,y())(9-12)
ynhK*
**
其中K叫做y(x)在[Xn,Xn1]上的平均斜率。
对于平均斜率K,只要提供一种计算方法,公式(9-12)就给出一种数值解公式。
例如,用K1f(Xn,yn)计算K*,就得到一阶精度的欧拉公式;用K2f(Xn1,yn1)替代K*,就得到隐式欧拉公式;如果用K1,K2的算术平均值计算K*,则可得到二阶精度的梯形公式。
由
此可以设想,如果在[Xn,xn1]上能多预测几个点的斜率值,用它们的加权平均值来计算K*,就有望得到具
有较高精度的数值公式,这就是Runge-Kutta法的基本思想。
Runge-Kutta公式的一般形式是
s
Kif(Xncih,ynhaijKj)
s'(i1,L,s)(9-13)
yn1ynhbiKi
i1
其中h为步长,s称为方法的级数,Ki是yy(x)在xnqh(0ci1)点的斜率预测值,G,b,aj均为常数。
这些常数的选取原则是使公式(9-13)具
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- v2 微分方程 数值 wb