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）具
升级会员 