第十五章常微分方程的解法文档格式.docx
《第十五章常微分方程的解法文档格式.docx》由会员分享，可在线阅读，更多相关《第十五章常微分方程的解法文档格式.docx（24页珍藏版）》请在冰豆网上搜索。
代替
代入
（1）中的微分方程，则得
化简得
如果用
的近似值
代入上式右端，所得结果作为
的近似值，记为
，则有
（2）
这样，问题
（1）的近似解可通过求解下述问题
（3）
得到，按式（3）由初值
可逐次算出
式（3）是个离散化的问题，称为差分方程初值问题。
需要说明的是，用不同的差商近似导数，将得到不同的计算公式。
）用数值积分方法
将问题
（1）的解表成积分形式，用数值积分方法离散化。
例如，对微分方程两端积分，得
（4）
右边的积分用矩形公式或梯形公式计算。
）Taylor多项式近似
将函数
在
处展开，取一次Taylor多项式近似，则得
再将
，得到离散化的计算公式
以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。
其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。
2欧拉（Euler）方法
2.1Euler方法
Euler方法就是用差分方程初值问题（3）的解来近似微分方程初值问题
（1）的解，即由公式（3）依次算出
这组公式求问题
（1）的数值解称为向前Euler公式。
如果在微分方程离散化时，用向后差商代替导数，即
，则得计算公式
（5）
用这组公式求问题
（1）的数值解称为向后Euler公式。
向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。
向前Euler公式是显式的，可直接求解。
向后Euler公式的右端含有
，因此是隐式公式，一般要用迭代法求解，迭代公式通常为
（6）
2.2Euler方法的误差估计
对于向前Euler公式（3）我们看到，当
时公式右端的
都是近似的，所以用它计算的
会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。
假定用（3）式时右端的
没有误差，即
，那么由此算出
（7）
局部截断误差指的是，按（7）式计算由
这一步的计算值
与精确值
之差
为了估计它，由Taylor展开得到的精确值
是
（8）
（7）、（8）两式相减（注意到
）得
（9）
即局部截断误差是
阶的，而数值算法的精度定义为：
若一种算法的局部截断误差为
，则称该算法具有
阶精度。
显然
越大，方法的精度越高。
式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。
3改进的Euler方法
3.1梯形公式
利用数值积分方法将微分方程离散化时，若用梯形公式计算式（4）中之右端积分，即
并用
这就是求解初值问题
（1）的梯形公式。
直观上容易看出，用梯形公式计算数值积分要比矩形公式好。
梯形公式为二阶方法。
梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为
（10）
由于函数
关于
满足Lipschitz条件，容易看出
其中
为Lipschitz常数。
因此，当
时，迭代收敛。
但这样做计算量较大。
如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler法。
3.2改进Euler法
按式（5）计算问题
（1）的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：
先用Euler公式求
的一个初步近似值
，称为预测值，然后用梯形公式校正求得近似值
，即
（11）
式（11）称为由Euler公式和梯形公式得到的预测—校正系统，也叫改进Euler法。
为便于编制程序上机，式（11）常改写成
（12）
改进Euler法是二阶方法。
4龙格—库塔（Runge—Kutta）方法
回到Euler方法的基本思想—用差商代替导数—上来。
实际上，按照微分中值定理应有
注意到方程
就有
（13）
不妨记
，称为区间
上的平均斜率。
可见给出一种斜率
，（13）式就对应地导出一种算法。
向前Euler公式简单地取
为
，精度自然很低。
改进的Euler公式可理解为
取
，
的平均值，其中
，这种处理提高了精度。
如上分析启示我们，在区间
内多取几个点，将它们的斜率加权平均作为
，就有可能构造出精度更高的计算公式。
这就是龙格—库塔方法的基本思想。
4.1首先不妨在区间
内仍取2个点，仿照（13）式用以下形式试一下
（14）
为待定系数，看看如何确定它们使（14）式的精度尽量高。
为此我们分析局部截断误差
，因为
，所以（14）可以化为
（15）
在点
作了Taylor展开。
（15）式又可表为
注意到
中
，可见为使误差
，只须令
（16）
待定系数满足（16）的（15）式称为2阶龙格—库塔公式。
由于（16）式有4个未知数而只有3个方程，所以解不唯一。
不难发现，若令
，
，即为改进的Euler公式。
可以证明，在
内只取2点的龙格—库塔公式精度最高为2阶。
4.24阶龙格—库塔公式
要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：
（17）
其中待定系数
共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差
的11个方程。
取既满足这些方程、又较简单的一组
，可得
（18）
这就是常用的4阶龙格—库塔方法（简称RK方法）。
5线性多步法
以上所介绍的各种数值解法都是单步法，这是因为它们在计算
时，都只用到前一步的值
，单步法的一般形式是
（19）
称为增量函数，例如Euler方法的增量函数为
，改进Euler法的增量函数为
如何通过较多地利用前面的已知信息，如
，来构造高精度的算法计算
，这就是多步法的基本思想。
经常使用的是线性多步法。
让我们试验一下
，即利用
计算
的情况。
从用数值积分方法离散化方程的（4）式
出发，记
，式中被积函数
用二节点
的插值公式得到（因
，所以是外插。
（20）
此式在区间
上积分可得
于是得到
（21）
注意到插值公式（20）的误差项含因子
，在区间
上积分后得出
，故公式（21）的局部截断误差为
，精度比向前Euler公式提高1阶。
若取
可以用类似的方法推导公式，如对于
有
（22）
其局部截断误差为
如果将上面代替被积函数
用的插值公式由外插改为内插，可进一步减小误差。
内插法用的是
，取
时得到的是梯形公式，取
时可得
（23）
与（22）式相比，虽然其局部截断误差仍为
，但因它的各项系数（绝对值）大为减小，误差还是小了。
当然，（23）式右端的
未知，需要如同向后Euler公式一样，用迭代或校正的办法处理。
6一阶微分方程组与高阶微分方程的数值解法
6.1一阶微分方程组的数值解法
设有一阶微分方程组的初值问题
（24）
若记
，则初值问题（24）可写成如下向量形式
（25）
如果向量函数
在区域
：
满足Lipschitz条件，即存在
，使得对
，都有
那么问题（25）在
上存在唯一解
问题（25）与
（1）形式上完全相同，故对初值问题
（1）所建立的各种数值解法可全部用于求解问题（25）。
6.2高阶微分方程的数值解法
高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。
设有
阶常微分方程初值问题
（26）
引入新变量
，问题（26）就化为一阶微分方程初值问题
（27）
然后用6.1中的数值方法求解问题（27），就可以得到问题（26）的数值解。
最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。
具体地说，对一阶线性微分方程组
（28）
阶方阵。
若矩阵
的特征值
满足关系
则称方程组（28）为刚性方程组或Stiff方程组，称数
为刚性比。
对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。
理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长
不作任何限制。
7Matlab解法
7.1Matlab数值解
7.1.1非刚性常微分方程的解法
Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。
Matlab的工具箱中没有Euler方法的功能函数。
）对简单的一阶方程的初值问题
改进的Euler公式为
我们自己编写改进的Euler方法函数eulerpro.m如下：
function[x,y]=eulerpro（fun,x0,xfinal,y0,n）;
ifnargin<
5,n=50;
end
h=（xfinal-x0）/n;
x
（1）=x0;
y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y1=y（i）+h*feval（fun,x（i）,y（i））;
y2=y（i）+h*feval（fun,x（i+1）,y1）;
y（i+1）=（y1+y2）/2;
例1用改进的Euler方法求解
解编写函数文件doty.m如下：
functionf=doty（x,y）;
f=-2*y+2*x^2+2*x;
在Matlab命令窗口输入：
[x,y]=eulerpro（'
doty'
0,0.5,1,10）
即可求得数值解。
）ode23，ode45，ode113的使用
Matlab的函数形式是
[t,y]=solver（'
F'
tspan，y0）
这里solver为ode45，ode23，ode113，输入
升级会员 