《微分方程数值解》由会员分享，可在线阅读，更多相关《微分方程数值解（12页珍藏版）》请在金锄头文库上搜索。
1、第四章 微分方程数值解4.1 算 法 当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程. 键入：syms x y %定义符号变量diff_equ= D2y+2*Dy-y=0; %D2y表示Dy=y=dsolve (diff_equ, x) %定义x为自变量y=cl*exp (2(1/2)-1)*x+c2*exp (-(2(1/2)+1)*x)%表达式中含c1与c2，表示通解.%初始条件为y (0)=0,(0)=1时，按如下方式调用y=dsolve (diff_equ, y (0)=0, Dy (0)=1, x)y=1/4*2(1/2)*exp (2(1
2、/2)-1)*x) 1/4*2(1/2)*exp (-(2(1/2)+1)*x)%画出函数y=y (x)的图形ezplot (y,-2,2) 图形具体形式请上机试之. 在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab求数值解的方法及应注意的问题. 例1 求解范德堡（vander pol）方程求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换则令 编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为待求解方程的函数名.m，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的文件. 文件存盘名为“vd
3、pol.m”.function yprime=vdpolyprime (1)=y (2); mu=2yprime (2)=mu*;yprime=yprime (1);yprime (2);说明 函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本例中，为解向量，为导数向量. yprime，yprime，函数返回vander pol方程的导数列向量. 因为所求结果为方程数值解，所以各向量维数只有在主程序求解时定下精度后才能确定. 主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外.clear functions%调试程序时，放置这一语句是必要的.
4、 它清除前边已编译的存在于内存中的废弃程序=ode23 (vdpol,0,30,1,0);y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( _ _)说明 龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为： =ode23 (,options) =ode45 (,options)其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. （1）若令，则输出在指定时刻给出，当时，输出在区间的等分点上给出，为步长. （2）若为自变量初值，为终值，此时，options决定自变量的维数，中的时间点不是等间隔的，这是为了保证所需的相对精度，
5、积分算法改变了步长. 用于设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若自行设定误差限，可用如下语句： options=odeset (reltol, abstol,)这里的与分别为设定的相对与绝对误差. 须注意的是无论用哪种方法确定的取值方式，必须由使用者确定且应与相匹配. 为初始条件，本例中，因为，这意味着解曲线，一般说，当解个未知函数的方程组时，为维向量，共含有个初始条件. 两个输出参数是列向量与矩阵，它们具有相同的行数，而矩阵的列数等于方程组的个数，本例中的列数为2，其中，为自变量上各点函数值，为上各点导数值. 最后，提请读者注意的是：ode45也不总是比
6、ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件 （4.1）其中适当光滑，关于满足Lipschitz条件，即存在使则（4.1）式的解存在且惟一. 关于的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点上寻求其数值近似解. 相邻两个节点间的间距称为步长，一般地取等步长，则 1、欧拉方法 在区间上用差商代替（4.1）式中，对中在上取值还是，而形成向前欧拉公式与
7、向后欧拉公式. （1）向前欧拉公式 取左端点，得如下公式 （4.2） 从点出发，由初值代入（4.2）求得 （4.3）反复利用（4.2），有 （4.4）误差其几何意义如图4.1所示. 图中为方程（4.1）的精确解曲线，其上任意点处切线斜率为. 从初值点出发，用该点斜率作一直线段，在处得到，由（4.2）式确定，再从出发，以为斜率作直线段，在处得到，这一过程继续下去，形成折线，作为积分曲线的近似，用 图4.1表示在处的精确值，为解的近似值，不难得到这一误差称为局部截断误差. 若一种算法局部截断误差为，则称该算法具有阶精度，所以向前欧拉公式具有1阶精度. （2）向后欧拉公式 若中取中的，则有如下公式：
8、 （4.5） 称式（4.5）为向后欧拉公式，因为此式中未知，故称其为隐式公式，无法用其直接计算，一般用向前欧拉公式产生初值.再按下式迭代其误差估计如下精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2. 为讨论局部截断误差，在图4.2中设点落在积分曲线上，按式（4.4）及式（4.5）分别得点为与，且点一定在积分曲线上相应点的上、下两边,所以将式（4.4）与（4.5）平均之，一定能得到更好的结果. （3）梯形公式 将向前与向后欧拉公式加以平均得到所 图4.2谓梯形公式 （4.6）其局部截断误差为，具有2阶精度. （4）改进的欧拉
9、公式 为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步 （4.7）或写为 （4.8） 最后指出，上述欧拉方法可推广至微分方程组，如向前欧拉公式为 2、龙格_库塔方法 由微分中值定理又因为，所以从而有 （4.9）令，称其为区间上的平均斜率，由（4.9）可知，给出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧拉公式中取，精度提高，下面，我们在区间内多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式 （4.10）其中，由于4个未知数只有3个方程，所以解不惟一，若令，即得改进的欧拉公式，具有2阶
10、精度. （3）4阶龙格_库塔公式 只给出精典格式中最常用的一种. （4.11）其计算精度为4阶4.3 模型与实验 1、模型与问题 例2 单摆运动 图4.3中一根长的细线，一端固定，另一端悬挂质量为的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 为平衡位置，在小球摆动过程中，当与平衡位置夹角为时，小球所受重力在其动运轨迹的分量为 （负号表示力的方向使减少），由牛顿第二定律可得微分方程 （4.12） 设小球初始偏离角度为，且初速为0，式（4.12）的初始条件为 （4.13） 当不大时，式
11、（4.12）化为线性常系数微分方程 图4.3 （4.14）解得 （4.15）简谐运动的周期为. 现在的问题是：当较大时，仍用近似，误差太大，式（4.12）又无解析解，试用数值方法在两种情况下求解，画出的图形，与近似解（4.15）比较，这里设. 例3 捕食与被捕食 当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，当甲独立生存时它的（相对）增长率与种群数量成正比，即有，为增长率，而乙的存在使甲的增长率减少，设减少率与乙的数量成正比，而得微分方程 （4.16）比例系数反映捕食者掠取食饵的能力. 乙离开甲无法生存，设乙独自存在时死亡率为，甲为乙提供食物，使乙的死亡率降低，而促其数量
12、增长，这一作用与甲的数量成正比，于是满足 （4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 （4.18）则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量、随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题： （1）设，求方程（4.16），（4.17）在条件（4.18）下的数值解，画出的图形及相图，观察解的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. （2）从式（4.16）和（4.17）消去得到 （4.19）解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数. （3）将方程（4.17）改写为 （4.20）在一个周期内积分，得到一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较. 2、实验 （1）方程（单摆问题）无解析解，为求其数值解，先将其化为等价的一阶方程组. 令，方程化为其中，为（弧度）与（弧度）两种情况，具体编程如下：先以danbai.m为文件名存放待解方程. 键入：function xdot =danbai (t,x) %x=x (1),x (2),g=9.8;1=25; %x (1)为解向量
