微分方程数值解.docx
《微分方程数值解.docx》由会员分享，可在线阅读，更多相关《微分方程数值解.docx（15页珍藏版）》请在冰豆网上搜索。
微分方程数值解
微分方程数值解及其应用
绪论
自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中，常常需要求解微分方程.但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下，只能用近似法求解，数值解法是一类重要的近似方法.本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这
些算法在处理来自生活实际问题中的应用，并结合MATLA软件，动手编程予以解决.
1微分方程的初值问题⑴
1.1预备知识
在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题
（1）
dyf（x,y）
dx
y（xo）yo
这里fx,y是矩形区域R:
xxoa,yy。
b上的连续函数.
对初值问题
（1）需要考虑以下问题：
方程是否一定有解呢？
若有解，有多少个解呢？
下面给出相关的概念与定理.
定义1Lipschitz条件[1][2]:
矩形区域R:
XX。
a,yy°|b上的连续函数
fx,y若满足：
存在常数L0，使得不等式fx,%fx,y2L%y?
对所有
x,y1,x,y2R都成立，则称fx,y在R上关于y满足Lipschitz条件•
定理1解的存在唯一性定理⑴⑶：
设f在区域Dx,yaxb,yR上连
续，关于y满足Lipschitz条件，则对任意的xoa,b,yoR，常微分方程初值问题
（1）
当xa,b时存在唯一的连续解yx
该定理保证若一个函数fx,y关于y满足Lipschitz条件，它所对应的微分方程的初值问题就有唯一解•在解的存在唯一性得到保证的前提下，自然要考虑方程的求
解问题.求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法.
定义2微分方程数值解：
对初值问题
（1）寻求数值解就是寻求解yx在一
x1x2L
Xn
Xn1L
系列离散节点
上的近似解y°,yi,y2丄，yn,yni,L，相邻两个节点的间距hXn1Xn称为步长.在一般情况下假定hhi0,1,L为常数，这时节点为xnx0nh,n0,1,2丄.
要求微分方程数值解，首先要建立数值算法，即对初值问题
（1）中的方程离散
化，建立求解数值解法的递推公式.一类是计算ym时只用到前一点的值yn，称为单
步法；另一类是用到yn1前面k点的值yn,yn1,L,ynk1称为k步法.
对初值问题
（1）式的单步法可用一般形式表示为
yn1ynh（Xn,yn,yn1,h），
其中多元函数与fx,y有关，当含有yn1时，方法是隐式的；若中不含yn1，
则为显式方法，所以显式单步法可表示为
yn1ynh（Xn,yn,h）.
（2）
设yx是初值问题
（1）的准确解，称Tn1yxn1yxnh（xn,yxn,h）为
显式单步法
（2）的局部截断误差.若存在最大正整数p,使显式单步法
（2）式的局部截断误差满足Tn1yxhyxhx,y,hOhp1，则称
（2）式有p阶精度.1.2几种常用的数值解法及其分析、比较
1.2.1欧拉法与后退欧拉法
1）欧拉法：
欧拉曾简单地用差分代替微分，即利用公式
y（Xn1）hy（Xn）y（Xn）f（Xn,Yn）
h
将初值问题
（1）离散化，则问题
（1）可化为
（3）
yn1ynhf（Xn,yn）,焉X。
nh,
此方法称为欧拉法.
欧拉方法的几何意义在数值计算公式中体现了出来
.在xy平面上，一阶微分方
程的解yyx称作它的积分曲线.积分曲线上一点
x,y的切线斜率等于函数
fx,y，按函数fx,y在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致•
基于上述几何解释，从初始点Po（Xo,yo）出发，先依方向场在该点的方向上推进到
xXi上一点R，再从R依方向场的方向推进到xX2上一点F2，循环前进便作出一
条折线RFP2L，因此欧拉方法又称为折线法.若初值yo已知，则由（3）式可逐步算
出
yiyohf（xo,yo）,
y2yihf（Xi,yJ,
M
为了分析计算公式的精确度，通常可用泰勒展开将yXni在Xn处展开，则有
容易看出，欧拉法（3）式具有一阶精度.
如果（4）式右端积分用右矩形公式hfxn1,yxn1近似，则得到另一个公式
（5）
yn1ynhfXn1,yn1，
称为后退欧拉法.
值得一提的是：
后退欧拉法与欧拉公式有着本质的区别，后者是关于yn1的直接
计算公式，它是显式的，而（5）式的右端含有关于yn1的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化•具体迭代过程如
下：
首先利用欧拉公式丫补铁ynhf（Xn,yn）给出迭代初值丫补铁，把它代入（5）式的右端，使之转化为显式，直接计算得y
（1）ynhf（xn1,y（0））.如此反复进行，得
ynki1）ynhf（Xni,ynk1）k0,1,L，
则得到后退欧拉法的迭代公式
yn0）1ynhf（Xn,yn）
（k1）（k）、
yn1ynhf（Xn1,yn1）
可以看出，后退欧拉法具有一阶精度，且计算比较麻烦•
1.2.2梯形方法
为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用yn代替yXn，yn1代替yXn1，则得
Yn1Yn2fXn,YnfXn1,Yn1，（6）
称其为梯形方法•
梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为yn0）1ynhf（Xn，yn）
（k1）h（k）.（7）
yn1yn-fXn,%fX.1,%1
为了分析梯形公式的收敛性，将（6）与（7）式相减，得
yn1ynk11）hfXn1,yn1f，k0,1,2,L
2
因为fX,y满足Lipschitz条件，于是有yn1丫二1吐yn1ynk1，其中L为fX,y
2
关于y的Lipschitz常数.如果选取h充分小，使得匹1，则当k时有
2
ynk11）yn1，这说明迭代过程（7）式是收敛的⑷.容易推导得出梯形法（7）式是二
阶方法.
经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的•从上述的迭代
公式可以看出，每迭代一次都要重新计算fx,y的值，而且迭代又要进行若干次，
计算相当的复杂•为此，有没有比较简便的计算方法呢？
下面给出改进的欧拉方法•
1.2.3改进的欧拉方法
由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解
yni，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下
h
yn1yn-fXn,ynfX.仆y.hf（8）
2
改进的欧拉法与梯形法一样，是二阶方法.
1.2.4Runge-Kutta方法
由前面讨论可知，从（4）式
Xn1
yXn1yXnXfX,yXdx
Xn
可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必
然要增加求积积点，为此将（4）式的右端用求积公式表示为
Xn1r
fx,yxdxhfXnih,yXnih，（9）
Xni1
一般来说，点数r越多，精度越高，上式右端相当于增量函数x,y,h，为得到便于
计算的显式方法，将公式（9）表示为：
yn1ynhXn,yn,h,（10）
其中
r
Xn,yn,hCjKi
i1
K1fXn,yn（11）
i1
KifXnih,ynhjKj,i2,Lr
j1
这里c,i,j均为常数.c为加权因子，Ki为第i段斜率，共有r段.我们把（10）和
（11）称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用
的一个公式：
K4fXnh,ynhK3
Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值
是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算ym，只用
到前面一步的值yn即可，因此每步的步长可以独立取定•常用的Runge-Kutta方法精度较
高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h可取得大些，求解区间
上的总步数可以少些.但Runge-Kutta方法也有些缺点，比如四阶Runge-Kutta方法每
算一步需要四次计算fx,y的值，计算量较大（对于复杂的fx,y而言）.
2数值方法的应用实例网
例1对于初值问题y10y，分别用欧拉法、改进的欧拉法，梯形法求y1的
y01
近似值.
解：
易得该方程的解析解yxe10x，y14.5400e-005,为比较，将按不同
数值计算方法所得结果列表如下：
表1三种不同方法的数值结果
h
欧拉法
改进的欧拉法
梯形法
0.2
-1
1
1
0.1
0
9.7656E-004
5.4994E-005
0.01
2.6561E-005
4.6223E-005
4.5026E-005
0.001
4.3717E-005
4.5408E-005
4.5396E-005
0.0001
4.5173E-005
4.5400E-005
4.5400E-005
0.02
0.015
0.01
0.005
0
欧拉法误差
改进欧拉法误差梯形法误差
图1三种不同方法数值解与精确解的误差曲线
从表1中可以看出：
当h0.2时，三种方法均不稳定，计算结果严重偏离精确值；h0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当h0.01时，三种方法均稳定，但精确度有区别•可以看出，h越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h值.
图1反映的步长h0.01时，三种数值方法的所得数值解与解析解在［0,1］区间的误差
曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法.
例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题
dy
dx
0,1
2x
y,Xy
并比较三种方法的结果.
解：
方程为n1的伯努利方程，可求得解析解为
现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果:
图2（a）步长为0.2时R-K法和解析解比较
图2（b）步长为0.2时改进的Euler法和解析解比较
图2（c）步长为0.2时欧拉法和解析解比较
上图2（a），（b），（c）描述的是步长为0.2时，用欧拉法、改进的欧拉法，Runge-Kutta
法求解方程所得的数值解与解析解之间的对比图•由图可知，Runge-Kutta法所得数
值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计
升级会员 