《实验4微分方程数值解.ppt》由会员分享，可在线阅读，更多相关《实验4微分方程数值解.ppt（33页珍藏版）》请在三一办公上搜索。
1、数学实验,常微分方程数值解,实验4,1 欧拉方法和龙格-库塔(Runge-Kutta)方法,设有常微分方程初值问题：,常微分方程数值解：计算（精确）解y(x)在一系列离散点 x0 x1 x2 xn 上的近似值，通常选取相等的计算步长 h。,常微分方程结论:f 适当光滑，对y 满足Lipschitz条件,则方程存在唯一解.,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
2、(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.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.73
3、21,1 欧拉方法,1.1 向前欧拉公式,设初值为 y(x0)=y0，,yn+1=yn+h f(xn,yn),n=0,1,2,这就叫向前欧拉公式。其精度为1 阶（局部截断误差为O(h2)。,p1,p2,p3,p4,function dy=ode121(x,y)dy=y-2*x/y;return,1.3 梯形公式,梯形公式:,kn=f(xn,yn),kn+1=f(xn+1,yn+1),Pn,Pn+1,Qn+1,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
4、.5000 1.4164 0.6000 1.4860 0.7000 1.5525 0.8000 1.6165 0.9000 1.6782 1.0000 1.7379,function 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)=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 改进的欧拉公式,称为改进的欧拉公式，它可
5、写作,例 用改进的欧拉公式解方程,function dy=ode121(x,y)dy=y-2*x/y;return,改进的欧拉公式可以从几何上加以解释如下：,以Pn,Q两点处的斜率之平均值为斜率，求得Pn+1,欧拉方法可以推广到解高阶方程和方程组，如：,向前的欧拉公式为：,解线性方程组,可采用向量函数表示.,对于高阶方程,需先降阶化为一阶常微分方程组,然后再按上方法求解.如对高阶线性方程,转换为一阶常微分线性方程组:,注意到,这里l1,l2为加权系数a(0 a1),3个数均为待定常数。,l1,l2，a 的选取标准：让精度更高,2 龙格-库塔方法,由二元函数的泰勒展式（推导在教材P72）,得,当
6、满足l1+l2=1,2a l2=1 时，精度达到2阶。显然解不唯一。,特别当l1=l2=1/2，a=1时就是改进的欧拉公式。,满足l1+l2=1,2a l2=1 的公式,（39）,称为二阶龙格-库塔公式。,可以证明，当在小区间 xn,xn+1内只取2 点的二 阶龙格-库塔公式的最高精度为2 阶（O(h3)）。,类似地，在小区间选四个点，适当选取加权系数，得到精度为四阶O(h5)的四阶龙格-库塔公式：,同样可推广到解微分方程组和高阶方程。,3 四阶龙格-库塔方法的MATLAB实现,先对所解方程或方程组（高阶方程）编写ODE文件（M函数文件，Ordinary Differential Equati
7、on）保存在MATLAB设置的搜索路径，然后在解方程的脚本文件中调用：,t,x=ode23(f,ts,x0,options),t,x=ode45(f,tspan,x0,options),或,options用于设置误差限（缺省时相对误差为10-3，绝对误差为10-6），用语句实现为：,options=odeset(reltol,rt,abstol,at),ts=a,b表示自变量的初值到终值的区间或ts=t0,t1,t2,tf,2 应用问题求解,例1、海上缉私 海防某部缉私艇上的雷达发现正东方向C海里处有一艘走私船正以一定速度向正北方向行驶，缉私艇立即以最大速度前往拦截。用雷达进行跟踪时，可保持缉
8、私艇的速度方向始终指向走私船。建立任意时刻缉私艇的位置和缉私艇航线的数学模型，确定缉私艇追上走私船的位置，求出追上的时间。,模型 设在t=0时刻缉私艇发现走私船，走私船以速度a平行于Y轴正向行驶，缉私艇按速度b以指向走私船的方向行驶(ba)，在任意时刻 t 缉私艇位于点P(x,y),则,即,先建立一个函数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,
9、ts,x0);t xplot(t,x);gridgtext(x(t);gtext(y(t);pausefigure;plot(x(:,1),x(:,2);gridgtext(x);gtext(y),再编写一个脚本M文件,M文件以ex4t1作为文件名存盘，保存在搜索路径中。,设 a=20,b=40,c=15;,设走私船在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.
10、0000 0.2000 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 缉私艇的位置 走私船的位置,例2、食饵-捕食者模型（弱肉强食）,问
11、题：自然界中不同物种之间存在着一种非常有趣的相互依存、相互制约的生存方式，种群甲靠丰富的自然资源生长，而种群乙靠捕食甲为生，鱼和鲨鱼、兔和山猫、落叶松和蚜虫都是这种生存方式的典型。,先作一些简化问题的假设，建立数学模型：,生态学上称种群甲为食饵(prey)，种群乙为捕食者(predator)，二者共处组成食铒捕食者系统（简称PP系统）。几十年来许多数学家和生态学家对这一系统进行了深入的研究，建立了一系列数学模型.下面介绍由意大利数学家Volterra提出的最简单的一个模型。,-d,甲为乙提供食物，使死亡率降低，并使其增长，设这个作用与甲的数量成正比，于是：,比例系数 a 反映捕食者掠取食饵的能
12、力。,捕食者独立存在时的死亡率为d，,比例系数b 反映食饵对捕食者的供养能力。,设食饵和捕食者的初始数量分别为：,x(0)=x0，y(0)=y0,食饵甲和捕食者乙在时刻 t 的数量(在一定区域的数量,密度)分别记作x(t),y(t)，当甲独立生存时，它的相对增长率为 r，,(5),(6),(7),r,x,而乙的存在使甲的增长率减小，设减小的程度与乙的数量成正比，于是 x(t)满足方程：,(,-a y),=r x a xy,y,(,+b x),=-d y+b xy,ts=0:0.1:15;x0=25,2;t,x=ode45(shier,ts,x0)result=t,xplot(t,x),grid
13、,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=diag(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文件定
14、义微分方程组。,用p_p.m调用ode45函数解方程组初值问题,绘制积分曲线和相图,食饵-捕食者模型的数值解,食饵-捕食者模型的相图,maxval=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,求最大、最小值,求周期与平均值,从前面所得的方程组初值问题,通过改写得到,(9),两边在一个周期内积分并除以周期T 得：,x(
15、t)和 y(t)的平均值刚好与相轨线的中心点P 的坐标相符。,上机内容,1、P84 1（1）（2）（3）龙格库塔方法2、P84 2，6,4 火箭发射,小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直发射时，燃料燃烧率为18kg/s，由此产生推力32000N，火箭引擎在燃料用尽时关闭。设火箭上升时，空气阻力与速度的平方成正比，比例系数为0.4kg/m。求引擎关闭瞬间火箭的高度、速度、加速度及火箭到达最高点时的高度和加速度，并画出高度、速度和加速度随时间变化的图形。,建立数学模型应考虑分有燃料和燃料已用尽两个阶段。设燃料燃烧率为常量，燃料燃烧的时间可以确定第一阶段的时间。第二阶段
16、的初值由第一阶段的终值确定。引擎关闭的时刻即第一阶段终止第二阶段开始的时刻。火箭先加速，后减速，达到最高点的时刻火箭的运动有什么特征？,第一阶段，燃料燃烧推动阶段，为变加速运动，初值问题为,M0=1400kg，Q0=1080kg，r=18kg/s，k=0.4kg/m，N=32000N,第二阶段，燃料烧尽，火箭作减速运动，初值问题为,求x=T 时，hT=x(T),求当,时的高度hm和加速度a。,建立数学模型，明确所要求出的量,function dx=odehj1(t,x)global M r N kg=9.8;dx=zeros(2,1);dx(1)=x(2);dx(2)=(N-(M-r*t)*g
17、-k*x(2)2)/(M-r*t);,function dx=odehj2(t,x)global M Q kg=9.8;dx=zeros(2,1);dx(1)=x(2);dx(2)=-g-k*x(2)2/(M-Q);,引擎关闭前的阶段：,引擎关闭后的阶段,clearglobal M Q r N k%声明全局变量M=1400;Q=1080;N=32000;r=18;k=0.4;T=Q/r;x0=0 0;t1,x=ode45(odehj1,0:T,x0);%求第一阶段解h1=x(end,1);v1=x(end,2);%引擎关闭时火箭高度和速度y0=h1,v1;%第二阶段初值t2,y=ode45(o
18、dehj2,T:T+20,y0);%求第二阶段解s=max(find(y(:,2)=0);%速度为零时的时间分量hm=y(s,1);%最大高度t=t1;t2(2:end);h=x(:,1);y(2:end,1);%高度与时间v=x(:,2);y(2:end,2);%速度函数,huojian.m,L=length(h);dt=(T+12)/(L-1);%用三点公式求加速度ac(1)=(-3*v(1)+4*v(2)-v(3)/2/dt;%起始点加速度for k=2:L-1 ac(k)=(v(k+1)-v(k-1)/2/dt;%中间点加速度end ac(L)=(v(L-2)-4*v(L-1)+3*v
19、(L)/2/dt;%最后点加速度a1=ac(length(t1)-1);%关闭引擎时加速度am=ac(length(t1)+s-1);%最高点加速度result1=T h1 v1 a1%引擎关闭时高度,速度和加速度result2=T+t2(s)hm am%达到最高点时间高度和加速度subplot(1,3,1),plot(t,h,b-),axis(0 80 0 14000)%高度函数subplot(1,3,2),plot(t,v,r-),axis(0 80-10 300)%速度函数subplot(1,3,3),plot(t,ac,g-),axis(0 80-30 40)%加速度函数,result
20、1=1.0e+004*0.0060 1.2190 0.0267 0.0002result2=1.0e+004*0.0071 1.3115-0.0005,实验内容：,P120 4.2 1)，3)，5)，其中3)题的建模已有提示，5)题如下图，建立坐标系后列出小船运动方程。注意不同情况下，小船到达对岸的时间要通过反复实验确定。,2 底端有孔容器的水面高度（实验内容之一）,问水从小孔中流完需要多少时间？2分钟时水面高度是多少？,r,先作适当简化假设，建立微分方程，求出其解，给出解的图形，回答上面的问题。,p r2dx=vp d02/4dt,考虑时间t 为水位x 的函数，作图再反过来。,x 0 0.1
21、 0.2 0.3,0.9 1.0 1.1 1.2,0.4 0.5 0.6 0.7 0.8,D 0.03 0.05 0.08 0.14,0.19 0.33 0.45 0.68 0.98,1.10 1.20 1.13 1.00,h,(单位m),h=1.2m d0=3cm v=0.6,容器直径D与高度x的关系如下表：,求水流完的时间和2分钟时的水面高度。,用三次样条插值求出容器边缘曲线函数确定水位x 处容器半径。,同学们再见！,同学们再见！,实验5 实验指导书 实验目的：1)掌握用MATLAB软件求微分方程数值解的欧拉方法和龙格-库塔方法;理解这两种方法的本质，除直接引用内部函数外，还能自编这两种方
22、法求解的M文件。2)掌握初步建模的知识，建立好几个简化的实际问题的微分方程数学模型。3)学习用微分方程数学模型解决简化的实际问题，通过对结果的分析，解释实际问题中的某些现象，预见事物发展的趋势，结果。实验内容：P120 4.2 1)a)c)，2)，3)，4)，5)。,2)对课堂讲授的例子，通过初始条件或参数的改变，产生了新的解答，从问题的结果比较中分析初始条件与参数对解的周期、食饵与捕食者的最大值、最小值、平均值的影响，作出在实际问题中的解释。注意参数的改变要在基本符合实际意义的范围内，要得到比较有说服力的结论，可能要多做几种实验比较。,指导意见：1)编写好欧拉公式与改进的欧拉公式求微分方程数
23、值解的程序M文件备用。对所给定的几个方程，给出了精确解，要求将你用两种方法计算的结果与精确解作出比较。c)为二阶方程，先化为等价的方程组初值问题，写出方程组的ODE函数文件，再用两种方法求解。对结果的比较分析作出评判。,5)*有两种建立坐标系的方法（以A为原点，y轴指向B和以B为原点，y轴指向A），建立的方程组有差异，初始条件也会有所不同。为避免在解方程时出现0做除数的情况，在做除数的量后加上eps这个特殊的常量。要得到小船到达彼岸的时间，需要反复实验。（本题的解析解存在，但比较复杂，不要求做）。,3)先作适当简化假设，以水位x为自变量，时间t 为未知函数的因变量建立微分方程，求出其解，给出解
24、的图形可将自变量因变量反过来。对第二问，须先用三次样条插值求得容器边缘曲线的函数，确定水位x处容器的半径。以下就和第一问一样建立微分方程模型求解了。微分方程以水位为自变量,时间为未知函数因变量。,一只小船从宽为d的河流的某岸A点正对着彼岸的B点渡河，已知河水流速v1 与船在静水中的速度v2 之比为k.a)建立小船航线的方程；b)设d为100m，v1=1 m/s，v2=2 m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图。c)求流速v1为0,0.5,1.5,2 m/s，v2不变，结果将如何？,A,B,v,v1,v2,水流方向,function dotx=odeguohe(t
25、,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),
