资料试验微分方程数值解.pptx

《资料试验微分方程数值解.pptx》由会员分享，可在线阅读，更多相关《资料试验微分方程数值解.pptx（32页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、例 用向前欧拉公式解方程function,z=foeula(f,a,b,y0,h)m=(b-a)/h;x(1)=a;y(1)=y0;for,n=1:m,x(n+1)=x(1)+n*h;,y(n+1)=y(n)+h*feval(f,x(n),y(n);endz=x,y;z=foeula(ode121,0,1,1,0.1);y1=sqrt(1+2*x);,sol,=z,y1sol=,0,1.0000,1.0000,0.1000,1.1000,1.0954,0.2000,1.1918,1.1832,0.3000,1.2774,1.2649,0.4000,1.3582,1.3416,0.5000,1.

2、4351,1.4142,0.6000,1.5090,1.4832,0.7000,1.5803,1.5492,0.8000,1.6498,1.6125,0.9000,1.7178,1.6733,1.0000,1.7848,1.73211,欧拉方法1.1,向前欧拉公式设初值为,y(x0),=,y0，,yn+1,=,yn,+,h,f,(,xn,yn,),n,=,0,1,2,这就叫向前欧拉公式。其精度为1,阶（局部截断误差为O,(h2)。,p1,p2,p3,p4x yox1x2x3x4function,dy=ode121(x,y)dy=y-2*x/y;return畴扎奥汉牺男呛蹄演垃坛珊懦蝴疲诚亨乍盈

3、诚涨探赶氏纵蓑取腕盾另臣替实验4微分方程数值解实验4微分方程数值解第1页/共32页y=adveula(ode121,0,1,1,10)y,=,0,1.0000,0.1000,1.0959,0.2000,1.1841,0.3000,1.2662,0.4000,1.3434,0.5000,1.4164,0.6000,1.4860,0.7000,1.5525,0.8000,1.6165,0.9000,1.6782,1.0000,1.7379function,s=adveula(f,a,b,y0,m)x=zeros(1,m+1);y=zeros(1,m+1);h=(b-a)/m;,x(1)=a;y(1

4、)=y0;for,n=1:m,x(n+1)=x(n)+h;,k1=feval(f,x(n),y(n);,k2=feval(f,x(n+1),y(n)+h*k1);,y(n+1)=y(n)+h*(k1+k2)/2;ends=x,y;1.4,改进的欧拉公式先由向前欧拉公式算出yn+1的预测值，再把它代替梯形公式的右端的yn+1，作为校正，即称为改进的欧拉公式，它可写作例,用改进的欧拉公式解方程function,dy=ode121(x,y)dy=y-2*x/y;return弃认蕴孰纠慕敷井胞权奢指辗姐刷缨却旨星倔椎讫贪砍寞峻栋谨茁七碟绰实验4微分方程数值解实验4微分方程数值解第3页/共32页改进的欧

5、拉公式可以从几何上加以解释如下：Q(xn+1,yn+hk1)k1k2 yxn xn+1Pn(xn,yn)y=f(x)xoPn+1(xn+1,yn+1)以Pn,Q两点处的斜率之平均值为斜率，求得Pn+1裕瑚捷厌器方每长盆烛降痈榨愿秃垒它勘规堵摇尖惺筐婆焕拣酷怠缀妥只实验4微分方程数值解实验4微分方程数值解第4页/共32页欧拉方法可以推广到解高阶方程和方程组，如：向前的欧拉公式为：解线性方程组,可采用向量函数表示.对于高阶方程,需先降阶化为一阶常微分方程组,然后再按上方法求解.如对高阶线性方程常用的方法是:令转换为一阶常微分线性方程组:阶曾贮灼逃晶末八钟捞痹寐勿瓶靳破射太专睁嘲打狐觉锋匣注株痴东类

6、宠实验4微分方程数值解实验4微分方程数值解第5页/共32页,注意到这里l1,l2为加权系数a,(0,aa)，在任意时刻,t,缉私艇位于点P(x,y),则即OXYc先建立一个函数M文件function,dx=jisi(t,x)a=20;b=40;c=15;s=sqrt(c-x(1)2+(a*t-x(2)2);dx=(c-x(1)/s;(a*t-x(2)/s*b;,%以列向量形式表示方程clearts=0:0.05:0.5;x0=0,0;t,x=ode45(jisi,ts,x0);t,xplot(t,x);gridgtext(x(t);gtext(y(t);pausefigure;plot(x(:

7、,1),x(:,2);gridgtext(x);gtext(y)再编写一个脚本M文件M文件以ex4t1作为文件名存盘，保存在搜索路径中。设,a=20,b=40,c=15;做离爵卖刃掘鳞畴眠蒂逗儡硒鄙枕娜陈咏普慕霍宋瑟咎妹如盼卿脐笔窑株实验4微分方程数值解实验4微分方程数值解第9页/共32页设走私船在t时刻的位置为x1(t),y1(t),则,0,0,0,15.0000,0,0.0500,1.9984,0.0698,15.0000,1.0000,0.1000,3.9854,0.2924,15.0000,2.0000,0.1500,5.9445,0.6906,15.0000,3.0000,0.200

8、0,7.8515,1.2899,15.0000,4.0000,0.2500,9.6705,2.1178,15.0000,5.0000,0.3000,11.3496,3.2005,15.0000,6.0000,0.3500,12.8170,4.5552,15.0000,7.0000,0.4000,13.9806,6.1773,15.0000,8.0000,0.4500,14.7451,8.0273,15.0000,9.0000,0.5000,15.0046,9.9979,15.0000,10.0000时刻t,缉私艇的位置,走私船的位置苞惜完锅殿己区宗侵终风屹伟伍抹练挣免蓟蹄滇儒犊第吓壹称歉巢贾疤

9、捆实验4微分方程数值解实验4微分方程数值解第10页/共32页,例2、食饵-捕食者模型（弱肉强食）问题：自然界中不同物种之间存在着一种非常有趣的相互依存、相互制约的生存方式，种群甲靠丰富的自然资源生长，而种群乙靠捕食甲为生，鱼和鲨鱼、兔和山猫、落叶松和蚜虫都是这种生存方式的典型。先作一些简化问题的假设，建立数学模型：,生态学上称种群甲为食饵(prey)，种群乙为捕食者(predator)，二者共处组成食铒捕食者系统（简称PP系统）。几十年来许多数学家和生态学家对这一系统进行了深入的研究，建立了一系列数学模型.下面介绍由意大利数学家Volterra提出的最简单的一个模型。租惑庄髓担涡节靶炸屎郭渴玩

10、饱肘校立泰硝秩畜篷译篡腿魂澡氓斤豆篡鲜实验4微分方程数值解实验4微分方程数值解第11页/共32页-d,甲为乙提供食物，使死亡率降低，并使其增长，设这个作用与甲的数量成正比，于是：比例系数,a,反映捕食者掠取食饵的能力。,捕食者独立存在时的死亡率为d，比例系数b,反映食饵对捕食者的供养能力。设食饵和捕食者的初始数量分别为：,x(0),=,x0，,y(0),=,y0,食饵甲和捕食者乙在时刻,t,的数量(在一定区域的数量,密度)分别记作x(t),y(t)，当甲独立生存时，它的相对增长率为,r,，(5)(6)(7)r,x,而乙的存在使甲的增长率减小，设减小的程度与乙的数量成正比，于是,x(t),满足方

11、程：(-a,y)=,r,x,a,xyy(+,b,x)=,-d,y,+,b,xy油蝎努质袍苫绎伯戈席越畸延湍推鳃弯玲瓜醚穗迈弯那掇韵匀磷佛箕孔谁实验4微分方程数值解实验4微分方程数值解第12页/共32页ts=0:0.1:15;x0=25,2;t,x=ode45(shier,ts,x0)result=t,xplot(t,x),grid,gtext(x1(t),gtext(x2(t),figure;plot(x(:,1),x(:,2),grid,.,xlabel(x1),ylabel(x2)function,xdot=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=dia

12、g(r-a*x(2),-d+b*x(1)*x;设r=1,d,=0.5,a=0.1,b=0.02,x0=25,y0=2，求初值问题的数值解。我们先将方程组改写成：初始条件改写为x(0),=,(,x10,x,20,)T,=,(,x0,y0,)T,上述方程组可写成：再令,x,=(,x1,x2,)T，方程组就成为下面的形式：用函数shier.m文件定义微分方程组。用p_p.m调用ode45函数解方程组初值问题绘制积分曲线和相图搬怒棵藻脊铆雅公熬俯毛笨或媚肝椽忱鹰纺标捅愤占挽炔专沮肾实岩佑念实验4微分方程数值解实验4微分方程数值解第13页/共32页食饵-捕食者模型的数值解食饵-捕食者模型的相图maxva

13、l=max(x)maxval,=,99.3910,28.4610minval=min(x)minval,=,2.0106,2.0000从图形和数值解中近似得周期T=10.7;x1p=trapz(x(1:108,1)*0.1/T;x2p=trapz(x(1:108,2)*0.1/T;x1p,x2pans,=,24.9776,10.0063求最大、最小值求周期与平均值怒稍鸯腕览伪轨该知掂村涟鄙遭铃冉涪垃耸沥匆紊坍佯绞氢阉孰宋萌婴拎实验4微分方程数值解实验4微分方程数值解第14页/共32页上机内容1、P84 1（1）（2）（3）龙格库塔方法2、P84 2，6卡结录郊腋裂诞潦杯赚腐兑级痊篆诉涸建兹缎仓

14、溜支躇务鹏眉讯枕镜说间实验4微分方程数值解实验4微分方程数值解第16页/共32页率它渊劫擎称除盆循酚鼠铀窑际雁亭党幕萌皇戌朗躇选鸯小沃俗哦大熙昔实验4微分方程数值解实验4微分方程数值解第26页/共32页实验5,实验指导书,实验目的：1),掌握用MATLAB软件求微分方程数值解的欧拉方法和龙格-库塔方法;理解这两种方法的本质，除直接引用内部函数外，还能自编这两种方法求解的M文件。2),掌握初步建模的知识，建立好几个简化的实际问题的微分方程数学模型。3),学习用微分方程数学模型解决简化的实际问题，通过对结果的分析，解释实际问题中的某些现象，预见事物发展的趋势，结果。实验内容：P120,4.2,1)

15、,a),c)，,2)，,3)，,4)，,5)。角烬萄杏京冰侥泳蛙阎尖堡洞推涪救峰击若聂尊享愚烙血充滓羡窘研先嚼实验4微分方程数值解实验4微分方程数值解第27页/共32页2),对课堂讲授的例子，通过初始条件或参数的改变，产生了新的解答，从问题的结果比较中分析初始条件与参数对解的周期、食饵与捕食者的最大值、最小值、平均值的影响，作出在实际问题中的解释。注意参数的改变要在基本符合实际意义的范围内，要得到比较有说服力的结论，可能要多做几种实验比较。指导意见：1),编写好欧拉公式与改进的欧拉公式求微分方程数值解的程序M文件备用。对所给定的几个方程，给出了精确解，要求将你用两种方法计算的结果与精确解作出比

16、较。c),为二阶方程，先化为等价的方程组初值问题，写出方程组的ODE函数文件，再用两种方法求解。对结果的比较分析作出评判。哇钨滁恬代侨缘断喜偏哲垒者啪并凑扑彻皂硫慷去置摊道特业扩斧噬琅蝎实验4微分方程数值解实验4微分方程数值解第28页/共32页5)*,有两种建立坐标系的方法（以A为原点，y轴指向B和以B为原点，y轴指向A），建立的方程组有差异，初始条件也会有所不同。为避免在解方程时出现0做除数的情况，在做除数的量后加上eps这个特殊的常量。要得到小船到达彼岸的时间，需要反复实验。（本题的解析解存在，但比较复杂，不要求做）。3),先作适当简化假设，以水位x为自变量，时间t,为未知函数的因变量建立

17、微分方程，求出其解，给出解的图形可将自变量因变量反过来。对第二问，须先用三次样条插值求得容器边缘曲线的函数，确定水位x处容器的半径。以下就和第一问一样建立微分方程模型求解了。微分方程以水位为自变量,时间为未知函数因变量。绚姨吊屈外题泵首峙跺恍肆炳搽盯痢料辐姆胺顽野娩诧哩盲抚辛窖廉鸯绣实验4微分方程数值解实验4微分方程数值解第29页/共32页,一只小船从宽为d的河流的某岸A点正对着彼岸的B点渡河，已知河水流速v1,与船在静水中的速度v2,之比为k,.,a),建立小船航线的方程；,b),设d为100m，v1=1,m/s，v2=2,m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作

18、图。,c),求流速v1为0,0.5,1.5,2,m/s，v2不变，结果将如何？ABvv1v2水流方向瑰部峪蹋捶侣涝蛆热加愧脆倪旺豌槐畅拳屉晚坦蔬画伸旱摄脊贝扰械疯磺实验4微分方程数值解实验4微分方程数值解第30页/共32页function,dotx=odeguohe,(t,x)global,v1,v2d=100;sq=sqrt(x(1)2+(d-x(2)2)+eps;dotx=zeros(2,1);dotx(1)=v1-v2*x(1)/sq;dotx(2)=v2*(d-x(2)/sq;根据数学模型,针对不同情况定义以下ODE函数，以备命令文件guohe.m调用。global,定义全局变量，由调用函数时输入，输入前也要用global申明所输入的为全局变量。命名为guohe.m的M文件,每次运行只要重新输入v1,v2的值即可得到不同结果到达彼岸的时间的确定要通过多次实验确定clearx0=0,0;global,v1,v2v1=1;v2=2;t,x=ode23(odeguohe,0,66.5,x0);t,x(:,1),x(:,2)plot(x(:,1),x(:,2)音珐团拟豺殴疲腊泼甲寄拥悦破材元嘴忠惕鳞有烦抚掀博赵聘蚁苞睹鞋舌实验4微分方程数值解实验4微分方程数值解第31页/共32页感谢您的观看！第32页/共32页