微分方程数值解.doc

第四章 微分方程数值解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)-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.首先编写待求解方程的文件. 文件存盘名为“vdpol.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%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序=ode23 (vdpol,0,30,1,0);y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( _ _)说明 龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为： =ode23 (,options) =ode45 (,options)其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. （1）若令，则输出在指定时刻给出，当时，输出在区间的等分点上给出，为步长. （2）若为自变量初值，为终值，此时，options决定自变量的维数，中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若自行设定误差限，可用如下语句： options=odeset (reltol, abstol,)这里的与分别为设定的相对与绝对误差. 须注意的是无论用哪种方法确定的取值方式，必须由使用者确定且应与相匹配. 为初始条件，本例中，因为，这意味着解曲线，一般说，当解个未知函数的方程组时，为维向量，共含有个初始条件. 两个输出参数是列向量与矩阵，它们具有相同的行数，而矩阵的列数等于方程组的个数，本例中的列数为2，其中，为自变量上各点函数值，为上各点导数值. 最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件 （4.1）其中适当光滑，关于满足Lipschitz条件，即存在使则（4.1）式的解存在且惟一. 关于的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点上寻求其数值近似解. 相邻两个节点间的间距称为步长，一般地取等步长，则 1、欧拉方法 在区间上用差商代替（4.1）式中，对中在上取值还是，而形成向前欧拉公式与向后欧拉公式. （1）向前欧拉公式 取左端点，得如下公式 （4.2） 从点出发，由初值代入（4.2）求得 （4.3）反复利用（4.2），有 （4.4）误差其几何意义如图4.1所示. 图中为方程（4.1）的精确解曲线，其上任意点处切线斜率为. 从初值点出发，用该点斜率作一直线段，在处得到，由（4.2）式确定，再从出发，以为斜率作直线段，在处得到，这一过程继续下去，形成折线，作为积分曲线的近似，用 图4.1表示在处的精确值，为解的近似值，不难得到这一误差称为局部截断误差. 若一种算法局部截断误差为，则称该算法具有阶精度，所以向前欧拉公式具有1阶精度. （2）向后欧拉公式 若中取中的，则有如下公式： （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）改进的欧拉公式 为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步 （4.7）或写为 （4.8） 最后指出，上述欧拉方法可推广至微分方程组，如向前欧拉公式为 2、龙格_库塔方法 由微分中值定理又因为，所以从而有 （4.9）令，称其为区间上的平均斜率，由（4.9）可知，给出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧拉公式中取，精度提高，下面，我们在区间内多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式 （4.10）其中，由于4个未知数只有3个方程，所以解不惟一，若令，即得改进的欧拉公式，具有2阶精度. （3）4阶龙格_库塔公式 只给出精典格式中最常用的一种. （4.11）其计算精度为4阶4.3 模型与实验 1、模型与问题 例2 单摆运动 图4.3中一根长的细线，一端固定，另一端悬挂质量为的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 为平衡位置，在小球摆动过程中，当与平衡位置夹角为时，小球所受重力在其动运轨迹的分量为 （负号表示力的方向使减少），由牛顿第二定律可得微分方程 （4.12） 设小球初始偏离角度为，且初速为0，式（4.12）的初始条件为 （4.13） 当不大时，式（4.12）化为线性常系数微分方程 图4.3 （4.14）解得 （4.15）简谐运动的周期为. 现在的问题是：当较大时，仍用近似，误差太大，式（4.12）又无解析解，试用数值方法在两种情况下求解，画出的图形，与近似解（4.15）比较，这里设. 例3 捕食与被捕食 当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，当甲独立生存时它的（相对）增长率与种群数量成正比，即有，为增长率，而乙的存在使甲的增长率减少，设减少率与乙的数量成正比，而得微分方程 （4.16）比例系数反映捕食者掠取食饵的能力. 乙离开甲无法生存，设乙独自存在时死亡率为，甲为乙提供食物，使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是满足 （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)为解向量, x (2)是导数xdot (1)=x (2); %xdot (1)=(1)=xdot (2)=-g/1*sin (x (1); %xdot (2)=xdot=xdot (1);xdot (2); %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下：clear functionst,x=ode23 (danbai,0,10,0.1745,0)%只计算（弧度）的情形.%对近似解，周期w=sqrt (9.8/25);y=0.1745*cos (w*t);t,x (:,1),y %显示数据，无分号.hold on %欲在同一幅图中画近似解.plot (t,x (:,1), b) %画数值解，绿色plot (t,y, g*) %用*号，红色画近似解.Hold off （2）食饵_捕食者 方程可化为如下形式初始条件表示为以shier.m存盘如下代码function xdot=shier (t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=diag (r-a*x (2),-d+b*x (1)*x;%x=,xdot=以main3.m存盘如下代码.clear functionsts=0:0.1:15;x0=25,2;t,x=ode45 (shier,ts,x0);t,x %显示数据t,x,yplot (t,x)gird %加网格线gtext (x(t),gtext (y(t), %用点鼠标方式pause,figure (2) %将x1(t),x2(t)放至指定点plot (x (:,1),x (:,2) %以x (1)与x (2)为坐标点画相图xlabel (x),ylabel (y) 可以猜测，与是周期函数，相图是封闭曲线，从数值解可近似定出周期为的最大和最小值分别为与的最大和最小值分别为和，为求与在一个周期的平均值，只需键入：y1=x (1:108,1); %x1周期为10.7.x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按y2=x (1:108,2); %梯形法对y1的积x2p=trapz (y2)*0.1/10.7, %分值可得 对方程化为积分得解 （4.21）即为原方程组的相轨线，其中由初始条件确定. 为说明上述相轨线是封闭的，令设其最大值分别为，若满足则有（令，可解出），又当时，相轨线（4.21）有意义. 当时，相轨线退化为一个点，对时，相轨线如图4.4（c），而图（a），（b）分别为与的图像，下面讨论如何由（a），（b）画（c）. （a） （b） （c）图4.4 设，若令，则有，由（a）知，使，且于是相轨线应过与，对，由，令，又由（b）知，存在使，且，于是相轨线又过与两点，所以对间每个，轨线总要通过纵坐标为的两点，于是相轨线是一条封闭曲线. 相轨线封闭等价于是周期函数，记周期为，为求其在一个周期内的平均值，由两边在一个周期内积分有：同理从而即的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与初始条件中的不是一件事. 注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比. 3、练习内容 （1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较. a）或. b) （Bessel方程，这里令，其精确解为）. c). （2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面积单位的小孔向外出水时，水柱