常微分方程的数值解法及其应用Word下载.doc
《常微分方程的数值解法及其应用Word下载.doc》由会员分享，可在线阅读，更多相关《常微分方程的数值解法及其应用Word下载.doc（49页珍藏版）》请在冰点文库上搜索。
1.1微分方程和解
含有未知量的等式称为方程，它表达了未知量所必须满足的某些条件。
方程是根据对未知量所进行的运算来分类的，如代数方程、超越方程等。
微分方程与代数方程和超越方程不同，它的未知量是函数，对其所施加的运算涉及求导或微分。
1.1.1微分方程的概念
一般说来，微分方程就是联系自变量、未知函数以及未知函数的某些导数或微分的关系式。
如果其中未知函数是一元函数，则称为常微分方程；
如果未知函数是多元函数，并且在方程中出现偏导数，则称为偏微分方程。
本论文主要介绍常微分方程，也简称微分方程或方程。
在一个微分方程中，出现的未知函数导数的最高阶数，称为方程的阶。
以为未知函数，为自变量的一阶常微分方程的一般形式可表示为
（1.1.1）
如果在（1.1.1）中能将解出，则得到方程
（1.1.2）
或
（1.1.3）
也称（1.1.1）为一阶隐式微分方程，（1.1.2）为一阶显式微分方程，（1.1.3）为一阶微分方程的微分形式。
阶隐式方程的一般形式为
（1.1.4）
阶显式方程的一般形式
（1.1.5）
方程（1.1.4）中，如果函数对未知函数和它的各阶导数都是一次的，则称其为线性常微分方程，否则，称其为非线性微分方程。
以为未知函数，为自变量的阶线性微分方程具有如下形式：
（1.1.6）
1.1.2微分方程的解——通解与特解
定义1.1设函数在区间上具有直到阶的导数。
如果把代入方程（1.1.4），有
在区间上关于恒成立，则称为方程（1.1.4）在区间上的一个解。
依据定义1.1可以直接验证：
（1）函数是方程在区间上的解，其中是任意常数。
另外，该方程还有两个解，它们不包含在前面解中。
（2）函数是方程在区间上的解，其中
和是独立的任意常数。
当然，都是方程的解，它们包含在前面解中。
从上面的讨论中看到事实：
微分方程的解可以包含任意常数，其任意常数的个数可以多到与方程的阶数相等，也可以不含任意常数。
阶常微分方程（1.1.4）的含有个独立的任意常数的解称为该方程的通解，而方程满足给定条件的解称为特解。
一般地，方程的特解可由其通解中任意常数取确定的常数导出。
以隐函数形式表示的通解称为通积分，而以隐函数形式表示的特解称为特积分，对于通解或者通积分的说法或使用，通常是不加区分的。
另外方程的通解不一定表示方程的所有解。
为了便于研究方程解的性质，常常需要考虑解的图像，或者以图形方式表示微分方程的解。
一阶微分方程的特解的函数图像是平面上的一条曲线，称为微分方程的积分曲线，而通解的函数图像是平面上的一族曲线，称为积分曲线族。
1.2常微分方程初值问题的一般提法
常微分方程初值问题的一般提法是求函数，满足
其中是已知函数，是已知值。
假设在区域上满足条件：
（1）在上连续；
（2）在上关于变量满足Lipschitz条件：
（1.2.3）
其中常数称为Lipschitz常数。
我们简称条件
（1）、
（2）的基本条件。
由常微分方程的基本理论，我们有：
定理1.1当在上满足基本条件时，一阶常微分方程初值问题（1.2.1）、（1.2.2）对任意给定存在唯一解在上连续可微。
定义1.2方程（1.2.1）、（1.2.2）的解称为适定的，若存在常数和，对任意满足条件及的和，常微分方程初值问题
（1.2.4）
存在唯一解，且。
适定问题的解连续依赖于（1.2.1）右端的和初值。
由常微分方程的基本理论，还有：
定理1.2当在上满足基本条件时，微分方程（1.2.1）、（1.2.2）的解是适定的。
在本章中假设在上满足基本条件，从而（1.2.1）、（1.2.2）的解存在且适定。
一般的一阶常微分方程组初值问题是求解
（1.2.5）
（1.2.5）的向量形式是
（1.2.6）
其中
记。
类似于定理1.1和定理1.2，有：
定理1.3若映射满足条件
（1）在上是从到上的连续映射；
（2）在上关于满足Lipschits条件；
任意。
则常微分方程组初值问题（1.2.5）存在的唯一的连续可微解而且解是适定的。
高阶常微分方程初值问题一般为
（1.2.7）
其中为给定值。
引进新的变量函数
（1.2.8）
则初值问题（1.2.7）化成了一阶常微分方程组初值问题
（1.2.9）
通过求解（1.2.9）得到（1.2.7）的解。
1.3初值问题数值解基本概念
初值问题的数值解法，是通过微分方程离散化而给出解在某些结点上的近似值。
在上引入结点
称为步长。
在多数情况下，采用等步长，即。
记（1.2.1），（1.2.2）的准确解为，记的近似值为，记为。
求值问题数值解的方法是逼近法，即在计算出后计算。
数值的方法有单步与多步法之分。
单步法在计算时只利用而多步法在计算时不仅要利用还有利用前面已算出的若干个。
我们称要用到的多步法为步法。
单步法可以看作多步法，但两者有很大差别。
步法只能用于的计算，要用其它的方法计算；
而且在稳定性上单步法比的多步法容易分析；
此外单步法容易改变步长。
单步法和多步法又都有显式方法和稳式方法之分。
单步显式法的计算公式可写成
（1.3.1）
隐式单步法的计算公式可写成
（1.3.2）
在（1.3.2）中右端项显含。
从而（1.3.2）是的方程式，要通过解方程求出。
显式多步法计算公式为
（1.3.3）
而隐式多步法计算公式为
（1.3.4）
右端项含。
多步法中一类常用方法是线性多步法
（1.3.5）
其中是独立于和的常数。
时（1.3.5）是显式的，时是隐式的。
数值解法及方法构造、误差分析等内容。
一些概念和定义在后面的论述中逐步引入。
第二章常微分方程的数值解法
一般说来通过初等积分法求解的方程仅有很少的一些类型，绝大部分从实际问题中提出的微分方程往往求不出其解析解。
而在实际问题中，对于复杂微分方程的求解，一般只需要得到解在若干个点上的近似值或者解的便于计算的近似表达式（只要满足所规定的精度）即可。
本章将选讲关于微分方初值问题的基本数值解法及其matlab程序设计。
2.1微分方程数值解法的一般步骤
一阶微分方程的初值问题为
（2.1）
寻求微分方程初值问题（2.1）的数值解，就是求解函数在一系列离散点（称为结点）上精确值的近似值。
在使用数值解法求解微分方程初值问题时，一般是按一下步骤：
（1）引入点列，其中，称为步长。
为了便于使用计算机进行编程计算，一般取步长为定值，即，
，（2.2）
（2）寻求数值解的方法，即寻求由计算出的递推公式。
（3）利用
（2）中的格式逐步求出近似解。
2.2Euler法
Euler法又称Euler折线法，它是解常微分方程的数值方法中最简单的一种方法。
Euler法的基本思想是：
在每一个小区间上，用一条切线来代替原函数曲线。
从整体上看就是用一条折线来代替解（或积分）曲线，并以此来求取一系列离散结点处函数近似值。
2.2.1Euler法的基本思想
对于一阶微分方程初值问题
考虑在结点的导数用差商代替，则（2.3）可近似地写成
，（2.5）
从出发，由初值（2.4），并利用（2.5）就可以得微分方程精确解的值的近似值，再以作为的近似值代入（2.5）的右端，可得的近似值，继续做下去，一般可写成
（2.6）
这种方法称为解初值问题的Euler法。
公式（2.6）称为求解初值问题（2.3）和（2.4）的Euler公式。
Euler法的几何意义如图3.1所示。
初值问题（2.3）和（2.4）的解曲线过点，从点出发以为斜率作直线，与直线交于点，显然，再从点出发以为斜率作直线段，与直线交于点，以此类推。
这样就得到解曲线的一条近似折线，所以，Euler法又称Euler折线法。
x
y
o
图2.1Euler法的几何意义
Euler折线法的递推结构：
（1）令，使用计算公式
计算；
（2）输出和，并使，重复以上过程。
2.2.2Euler法的局部截断误差进行估计。
在用Euler法求解一阶微分方程初值问题（2.3）和（2.4）时，实际上是利用Euler公式进行计算，计算出的近似值与精确值之间必然是有误差的，所产生的误差为
（2.7）
将函数在点处进行Taylor展开，
将代入上式，则
注意到
可得
（2.8）
观察（2.8），可以发现，当步长趋于零时，是的同阶无穷小，记为.由此可知，我们用Euler折线法求解一阶常微分方程初值问题（2.3）和（2.4）时，在某一结点处的局部截断误差与选择的步长的平方成正比。
如果求解公式的局部截断误差是，则称该求解公式是阶的，或者具有阶精度的，因此，Euler折线法是具有一阶精度的一种数值方法。
作为理论误差的应用，在用计算机编程进行数值解计算时，可以通过式（2.8）和给定的误差要求初步确定步长和选取的点个数。
例如，如果给定的计算误差为，则可通过
初步确定出步长，进行适当调整，使得等式成立，并尽可能地取为正数，为有限小数。
2.2.3Euler折线法的数值实例及matlab源程序
例2.1使用Euler折线法求解初值问题
取步长为。
解应用Euler折线法，在计算机中运用matlab软件，程序为：
clearclf
szy=[];
y=1;
szy=[szy;
y];
forx=0.1:
0.1:
1
y=y+（1+y）*0.1;
end
szy
输出为：
szy=1.00001.20001.42001.6620
1.92822.22102.54312.8974
3.28723.71594.1875
用折线法算出的初值问题在点处的数值解（存放在数组中），而输入：
y=dsolve（'
Dy=1+y'
'
y（0）=1'
）
可以算得此微分方程初值问题的精确解：
y=2*exp（x）–1
为了比较精确解和近似解的误差，编写程序算出精确解在点处的纵坐标（存放在数组jqy中）
jqy=[];
forx=0:
y=-1+2.*exp（x）;
jqy=[jqy;
jqy
现在将折线和曲线在同一坐标系中画出来，如图2.2所示．
x=0:
1;
plot（x,szy,'
o'
x,jqy）
图2.2折线法的数值解与函数真实值
在点处的精确解与近似解相除所得的差。
从输出容易观察到：
此方法在经过多步以后，其误差会积累起来．为了减少误差，一种方法是减少步长的值。
另外，在近似算法上也有很多改进办法。
2.2.4梯形法
可以注意到，在利用Euler法求解微分方程初值问题时计算简单、计算量小、程序简洁，但是，计算精度却很低。
为了提高计算精确度，同时又保持计算简单，程序简洁的优点，所以给出了梯形法。
其基本思想是：
在每个小区间上，仍然用一条直线来代替原函数曲线，但是这条直线较Euler折线法中的切线更近似于解函数曲线。
为寻求这条直线，考虑微分方程
在区间上对方程两边关于积分，则
于是有
即
，（2.9）
现在，利用左矩形公式来计算（2.9）中的积分，有
，
由此可得
.
如果用和分别代替和，则可看出这个迭代公式就是Euler公式。
如果改用梯形公式来计算（2.9）中的积分，即
代入（2.9）中，得到
在数值计算格式中，用和分别代替和，则得到
（2.10）
从出发，由初值,利用（2.10）可以依次得。
公式（2.10）是由数值积分中的梯形公式得出来的，将它称为梯形公式。
和前面的Euler公式相比。
梯形公式计算是也只用到了前一步的数值。
但是，如果已知，将代入（2.10）的右端，一般就不能直接得到，而需要通过其他方法（如迭代法）来求解，称梯形公式为单步隐式公式。
单步隐式公式的一般形式为
其中为增函数。
梯形公式的增量函数为
。
下面不讨论梯形法的局部截断误差。
将在处用Taylor展开，得:
将在处进行一阶Taylor展开，得:
将的一阶展开式代入梯形公式（2.10），得
因此，梯形法的误差为
由此可知，梯形公式的计算误差与步长的三次方称正比。
可见，梯形法具有二阶精度。
2.2.5改进的Euler法
比较梯形公式（2.10）和Euler公式（2.6），梯形公式是隐式方程，且所耗费的计算量大于Euler公式求出
用近似准确值，称为预测值。
然后用梯形公式将它校正为较为准确的值，称由此得到的为校正值，即
这样建立起来的预测-校正系统通常称为改进的Euler公式，即
（2.11）
为了便于编程计算，可将改进的Euler公式（2.11）改写为下列形式
（2.12）
改进的Euler法的递推结构为
（1）令，利用式（2.12）依次计算出和；
（2）输出，并使，重复以上过程。
2.2.6改进的Euler法求解微分初值问题实例及matlab源程序
例2.2用改进的Euler法求解例2.1中的初值问题，取。
解改进的此法可以算出此问题的数值解，存放在数组gjszy中。
程序为：
x=0:
h=0.1;
gjszy=[];
gjszy=[gjszy;
fori=1:
length（x）-1
y1=y+（1+y）*h;
y=y+（1+y+1+y1）*h/2;
end
gjszy
plot（x,gjszy,'
x,jqy,'
r'
比较结果如图2.3所示。
图2.3改进Euler方法的数值解与函数真实值
此方法求得的近似解与真实解几乎完全重合，与Euler法相比，改进的Euler法消除了误差，结果更为精确，就精确度而言，改进的Euler法比Euler法要好，它具有二阶精度。
2.3Runge-Kutta法
2.2节的Euler法和改进的Euler法的计算精度分别是一阶的和二阶的，现在的问题是：
能否构造更高阶精度的数值方法？
前面提到的几种数值解法的精度是很低的，下面给出高阶一步法——Runge-Kutta方法。
它是最常用的一种数值解法，因为它相当精确、稳定、容易编程。
Runge-Kutta方法至今仍然得到广泛地应用。
2.3.1.Runge-Kutta法的基本思想
Runge-Kutta法的一般形式是
，（2.13）
其中函数具有下列形式：
，（2.14）
，
（2.15）
其中，，等均为常数，为整数。
现在考虑方法的几何意义，Euler法是用斜率向前推进，而Runge-Kutta法则可以看成用个斜率的一种加权平均来向前推进，以达到提高计算方法的精度。
式（2.14）用到了个，称为阶的显式Runge-Kutta法。
对于差商，由微分中值定理，可得
由于，于是
（2.16）
称为区间上的平均斜率，记作。
因此，只要给出平均斜率的一种算法，再由式（2.16）就可以得到求解微分方程的一种数值计算公式。
对于Euler法和改进的Euler法，可发现在Euler公式（2.6）仅取一个点的斜率作为平均斜率的近似，计算精度较低。
而改进的Euler公式（2.11）则是利用了和两个点的斜率和的平均值作为平均斜率的近似，即，其中是通过已知的值利用Euler公式求得，此时计算精度提高到了二阶。
由上面的分析，找到了改进的Euler公式较Euler公式的计算更精确的原因。
那就是在确定平均斜率时，多取了一个点的斜率值。
因此，受此启发，如果设法在区间上多给出几个点的斜率值，然后加权平均得到的值作为的近似值，那么就有可能构造出精确度更高的数值计算公式，这就是构造Runge-Kutta法的基本思想。
2.3.2二阶Runge-Kutta法
构造二阶Runge-Kutta公式的主要步骤如下：
（1）在区间上取两点，给出相应的斜率值；
（2）对这两个斜率值加权平均作为平均斜率的近似值；
（3）对有关的函数进行Taylor展开，得到关于步长的幂函数。
为使计算公式达二阶，则其中的系数为零，从而建立了有关参数应满足的方程组；
（4）解此方程组就得到相应的二阶Runge-Kutta公式。
下面来详细说明构造二阶Runge-Kutta公式的主要步骤，在区间上取两点和，以这两点处的斜率值和加权平均来求取平均斜率的近似值，即
其中点的斜率值和点的斜率值分别为
。
根据假设，取
（2.17）
为计算公式。
其中和为待定常数。
对在处进行二阶Taylor展开，则有
（2.18）
对在处进行一阶Taylor展开，得
（2.19）
将式（2.18）和式（2.19）代入式（2.17）中，得
（2.20）
比较式（2.18）和式（2.20），得
（2.21）
方程组（2.21）中有三个未知量，却只有两个方程，因此它有无穷多组解。
如果去。
则有，用和分别代替和，此时计算公式（2.17）变为
这一数值方法就是改进的Euler公式。
凡是满足条件式（2.21）的计算格式为式（2.17）的这些格式统称为二阶Runge-Kutta格式。
改进的Euler公式就是众多二阶Runge-Kutta公式中的一种特殊格式。
现在给出其他二阶Runge-Kutta格式。
如果取，则，此时得到二阶Runge-Kutta法的计算公式为
（2.22）
其中点为区间的中点，因此计算公式（2.22）又称为变形的Euler公式。
2.3.3二阶Runge-Kutta法的简单数值实验
例2.3用二阶Runge-Kutta法求解求方程
解利用二阶Runge-Kutta公式（2.22）解给定的初值问题的形式为
MATLAB源程序为：
fun=inline（'
y^2+x^3'
x'
y'
）;
[x,y]=ode23（fun,x,0.
升级会员 