实验一 常微分方程

1、年级、专业 2012级数学与应用数学 姓名 张旭 学号 12160011063 名单序号 16 实验时间 2014年 3月 19 日 使用设备、软件 PC, MATLAB 注： 实验报告的最后一部分是实验小结与收获 实验一 常微分方程1. 分别用Euler法和ode45解下列常微分方程并与解析解比较： (1) 解析法：y=dsolve('Dy=x+y','y(0)=1','x')y = 2*exp(x) - x - 1Euler:function x,y=euler(odefun,xspan,y0,h)x=xspan(1):h:xspan(2);

2、y(1)=y0;for i=1:length(x)-1 y(k+1)=y(i)+h*feval(odefun,x(i),y(i); end x=x' ; y=y' ;endode45:odefun=inline('x+y','x','y');xspan=0,3;y0=1;h=0.1;x1,y1=euler(odefun,xspan,y0,h);x2,y2=ode45(odefun,xspan,y0);x3=0:0.1:3;y3=2*exp(x3)-x3-1;plot(x1,y1,'k',x2,y2,'ko&

3、#39;,x3,y3,'k*');xlabel('x轴');ylabel('y轴');legend('euler','ode45','dsolve');ode45求得的结果与用解析法求得的结果更接近，故ode45的精度较高，Euler法求得的结果精度较低。 (2) 令 则原方程等价于方程组：， ，不能解析，只能用数值法求解。Euler：function t,y=euler2(odefun1,odefun2,tspan,y0,h)t=tspan(1):h:tspan(2);y(1,1)=y0(1);y

4、(2,1)=y0(2);for i=1:length(t)-1 k1=odefun1(t(i),y(1,i),y(2,i); k2=odefun2(t(i),y(1,i),y(2,i); y(1,i+1)=y(1,i)+h*d1; y(2,i+1)=y(2,i)+h*d2; end t=t' y=y'endode45： odefun1=inline('0*t1+0*y1+y2'); odefun2=inline('-2*y1+0.01*y22+sin(t1)'); t1,y1=euler2(odefun1,odefun2,0,5,0,1,0.1)

5、; t2,y2=ode45('eu',0,5,0,1); plot(t1,y1(:,1),'o',t2,y2(:,1),'LineWidth',2); xlabel('t轴'); ylabel('y轴'); legend('euler','ode45');ode45中 eu：function dy=eu(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-2*y(1)+0.01*y(2)2+sin(t);ode45求得的结果精度较高，euler法求得的结果在准确值

6、上下波动。2. 一通过原点的曲线，它在处的切线斜率等于若上限增为1.58，1.60会发生什么？ 等价于求解 ，且的初值问题。解析法： y=dsolve('Dy=2*x+y2','y(0)=0','x')y =(2(1/3)*airy(3,-2(1/3)*x)+2(1/3)*3(1/2)*airy(1,-2(1/3)*x)/(3(1/2)*airy(0, -2(1/3)*x) + airy(2, -2(1/3)*x)ode45： odefun=inline('2*x+y2'); subplot(1,3,1);ode45(odefun

7、,0,1.57,0);title('0<x<1.57'); subplot(1,3,2);ode45(odefun,0,1.58,0);title('0<x<1.58'); subplot(1,3,3);ode45(odefun,0,1.60,0);title('0<x<1.60');曲线单调递增，且在x>1.5之后的斜率增长速度很快，若上限增为1.58，1.60，则相应的y将会出现更大的增长。3. 求解刚性方程组：function dy=fun(t,y)dy=zeros(2,1);dy(1)=-1000

8、.25*y(1)+999.75*y(2)+0.5;dy(2)=999.75*y(1)-1000.25*y(2)+0.5;ode45： t,y=ode15s('fun',0,50,1,-1); plot(t,y(:,1),'o',t,y(:,2),'k-','LineWidth',2); legend('y1','y2');4. (温度过程)夏天把开有空调的室内一支读数为20 的温度计放到户外，10分钟后读25.2 , 再过10分钟后读数28.32 。建立一个较合理的模型来推算户外温度。由题意可知由于

9、随着时间的增加，温度增长越来越慢，户外温度与温度计示数之差也越来越小，且温差为零时温度的增长率也为零，故可以认为温度的增长率与温差成正比，设户外温度为m，温度计的示数为y，比例系数为k，则可建立模型 解析法：y=dsolve('Dy=k*(c-y)','y(0)=20','t')y = m - (m - 20)/exp(k*t)由y（10）=25.2,y(20)=28.32建立方程组，消去k，得： （m-20）（m-28.32）=（m-25.2）（m-25.2）解得：m=33 所以户外温度约为33 。5. (广告效应)某公司生产一种耐用消费品，市

10、场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。(1) 建立该问题的数学模型，并求其数值解与模拟结果作以比较；设t时刻该消费品的市场占有率为y，建立方程：解析解：y=dsolve('Dy=0.5-0.5*y','y(0)=0.05')y = 1 - (19*exp(-t/2)/20数值解： odefun=inline('0.5-0.5*y','t','y'); t1,y1=ode45(odefun,0,10

11、,0.05); t2=0:0.1:10; y2=1-(19*exp(-t2/2)/20; plot(t1,y1,'o',t2,y2,'k'); legend('ode','dsolve');(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？由解析解可列出方程 1 - (19*exp(-t/2)/20=0.8，所以解得t=3.11636. (肿瘤生长) 肿瘤大小V生长的速率与V的a次方成正比，其中a为形状参数，0£a£1;而其比例系数K随时间减小，减小速率又与当时的K值成正比，比例系数为环境参数b。设某

12、肿瘤参数a=1, b=0.1, K的初始值为2，V的初始值为1。问(1) 此肿瘤生长不会超过多大？由已知列出方程组 ，代入具体数值，得 ，解析法：k,v=dsolve('Dk=-0.1*k','Dv=k*v','k(0)=2','v(0)=1','t')k =2*exp(-t/10)v =exp(20)*exp(-20*exp(-t/10)t=0:0.1:100; v=exp(20)*exp(-20*exp(-t/10); plot(t,v,'LineWidth',2); xlabel('t

13、轴'); ylabel('v轴');因肿瘤不断长大，故t趋于无穷时，该肿瘤达到最大，此时极限为exp(20)=4.8517*108，故此肿瘤生长不会超过4.8517*108 。(2) 过多长时间肿瘤大小翻一倍？令exp(20)*exp(-20*exp(-t/10)=2，解得t=-10*ln(1-1/20*ln2)=0.3527,(3) 何时肿瘤生长速率由递增转为递减？由已求得的结果可得与的关系为=2*exp(20-t/10)*exp(-20*exp(-t/10)，t1=0:0.1:100; v1=2*exp(20-t1/10).*exp(-20*exp(-t1/10);

14、 plot(t1,v1,'LineWidth',2); xlabel('t轴'); ylabel('v轴');显然，最大值处对应的t即为所求： t2=0:0.01:100; v2=2*exp(20-t1/10).*exp(-20*exp(-t1/10); m,n=max(v2); t=t2(n)得到t = 29.96 (4) 若参数a=2/3呢？1、建立方程组， 解析法：k,v=dsolve('Dk=-0.1*k','Dv=k*v(2/3)','k(0)=2','v(0)=1',&#

15、39;t')k = 2*exp(-t/10) 2*exp(-t/10) 2*exp(-t/10)v = -(20*exp(-t/10) - 23)3/27 (37/2 + (3(1/2)*3*i)/2 - 20*exp(-t/10)3/27 -(20*exp(-t/10) + (3(1/2)*3*i)/2 - 37/2)3/27取实解k = 2*exp(-t/10) v=-(20*exp(-t/10) - 23)3/27并画出v-t图像： t=0:0.1:100; v=-(20*exp(-t/10) - 23).3/27; plot(t,v,'LineWidth',2)

16、; xlabel('t轴'); ylabel('v轴');显然，当t趋于无穷时，该肿瘤达到最大，此时极限为-(- 23)3/27=450.6296，故此肿瘤生长不会超过450.6296 。2、令-(20*exp(-t/10) - 23)3/27=2，解得t=0.3977,3、由已求得的结果可得与的关系为=2*exp(-t/10)*(20*exp(-t/10) - 23)2/9，t1=0:0.1:100; v1=2*exp(-t1/10).*(20*exp(-t1/10) - 23).2/9; plot(t1,v1,'LineWidth',2);

17、xlabel('t轴'); ylabel('v轴');显然，最大值处对应的t即为所求，： t2=0:0.01:100; v2=2*exp(-t2/10).*(20*exp(-t2/10) - 23).2/9; m,n=max(v2); t=t2(n)t = 9.5900选做题：1.(生态系统的振荡现象)第一次世界大战中，因为战争很少捕鱼，按理战后应能捕到更多的鱼才是。可是大战后，在地中海却捕不到鲨鱼，因而渔民大惑不解。令x1为鱼饵的数量，x2为鲨鱼的数量，t为时间。微分方程为 式中a1, a2, b1, b2都是正常数。第一式鱼饵x1的增长速度大体上与x1成正比，即按a1x1比率增加, 而被鲨鱼吃掉的部分按b1x1x2的比率减少；第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a2x2的比率减少，但又根据鱼饵的量的变化按b2x1x2的比率增加。对a1=3, b1=2, a2=2.5, b2=1, x1(0)=x2(0)=1求解。画出解曲线图和相轨线图，可以观察到鱼饵和鲨鱼数量的周期振荡现象。代入具体数值后，原方程组变为： ，x1(0)=x2(0)=1 ode45