实验一常微分方程解析

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

2、);y(1)=y0;for i=1:le ngth(x)-1y(k+1)=y(i)+h*feval(odefun ,x(i),y(i);endx=x;y=y;endode45: odefu n=in li ne(x+y,x,y);xspa n=0,3;y0=1;h=0.1;x1,y1=euler(od efun,xspan,y0,h); x2,y2=od e45(od efun,xspan,y0); x3=0:0.1:3;y3=2*exp(x3)-x3-1;Plot(x1,y1,k,x2,y2,ko,x3,y3,k*); xlabel(x 轴);ylabel(y 轴);legend(euler

3、,od e45,dsolve);ode45求得的结果与用解析法求得的结果更接近，故Euler法求得的结果精度较低。3 / 92014春数学实验实验一常微分方程注:实验报告的最后一部分是实验小结与收获 y 0.01(y)2 2y =sin(t), y(0) = 0,y (0) = 1,0 : t 5令丫勺=y, y2 = y：则原方程等价于方程组：丿2y2 二 一2%0.01y2 sin(t)%(0) =0,y2(0) =1，0 :：: t 5，不能解析，只能用数值法求解。Euler: function t,y=eul er2(od efun 1,odefun2,tspan,y0,h) t=ts

4、pa n(1):h:tspa n( 2);y(1,1)=y0(1);y(2,1)=y0(2);for i=1:le ngth(t)-1k1=od efun 1(t(i),y(1,i),y(2,i);k2=od efun2(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;endt=t;y=y;endode45 :odefu n1=i nli ne(0*t1+0*y1+y2);odefu n2=i nlin e(-2*y1+0.01*y2A2+si n( t1); t1,y1=euler2(od efun 1,od efun

5、2,0,5,0,1,0.1); t2,y2=od e45(eu,0,5,0,1);plot(t1,y1(:,1),o,t2,y2(:,1),Li neWidth,2); xlabel(t 轴);ylabel(y 轴);legend(euler,od e45);ode45 中 eu: function dy=eu(t,y)实验时间2014年3月19日使用设备、软件PC, MATLAB注：实验报告的最后一部分是实验小结与收获dy=zeros(2,1);dy(1)=y(2);dy(2)=-2*y(1)+0.01*y(2)A2+s in (t);ode45求得的结果精度较高，euler法求得的结果在准

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

7、ot(1,3,1);od e45(od efun,0,1.57,0);title(0x1.57);subplot(1,3,2);od e45(od efun,0,1.58,0);title(0x1.58);subplot(1,3,3);od e45(od efun,0,1.60,0);title(0x1.60);Dix-tl .57口 a1.5之后的斜率增长速度很快，若x上限增为1.58 , 1.60，则相应的y将会出现更大的增长。年级、专业 2012级数学与应用数学姓名 张旭 学号 12160011063 名单序号16实验时间2014年3月19 日使用设备、软件PC, MATLAB实验一常微

8、分方程注:实验报告的最后一部分是实验小结与收获3. 求解刚性方程组:0 ： x ： 50.1000.25 y 999.75y20.5, %(0)=1,y2 =999.75力-1000.25 y20.5,y1(0 1,fun ctio n dy=fu n( t,y) dy=zeros(2,1); dy(1)=-1000.25*y(1)+999.75*y(2)+0.5; dy(2)=999.75*y(1)-1000.25*y(2)+0.5;ode45 :t,y=od e15s(fun,0,50,1,-1); plot(t,y(:,1),o,t,y(:,2),k-,Li neWidth,2); le

9、ge nd(y1,y2);5 / 92014春数学实验实验一常微分方程4. (温度过程)夏天把开有空调的室内一支读数为20 C的温度计放到户外，10分钟后读25.2 C,再过10分钟后读数28.32 C。建立一个较合理的模型来推算户外温度。由题意可知由于随着时间的增加，温度增长越来越慢，户外温度与温度计示数之差也越来越小，且温差为零时温度的增长率也为零，故可以认为温度的增长率与温差成正比，设户外温度为m，温度计的示数为 y，比例系数为k，则可建立模型y =k(m -y), y(0) =20解析法:y=dsolve(Dy=k*(c-y),y(0)=20,t)注:实验报告的最后一部分是实验小结与收

10、获y = m - (m - 20)/exp(k*t)由 y (10)=25.2,y(20)=28.32建立方程组m20“ cm 10k25-2em 20m -20k28-32e消去k，得:(m-20)(m-28.32)= ( m-25.2)( m-25.2 )解得：m=33所以户外温度约为 33C5. (广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后， 该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有 买的百分比成正比，且估得此比例系数为0.5。(1)建立该问题的数学模型，并求其数值解与模拟结果作以比较;设t时刻该消费品的市场占有率为y,

11、建立方程：/-0.5* (1 y), y(0) =5%解析解：y=dsolve(Dy=0.5-0.5*y,y(0)=0.05)y = 1 - (19*exp(-t/2)/20数值解：odefu n=inlin e(0.5-0.5*y,t,y);t1,y1=od e45(od efun,0,10,0.05); t2=0:0.1:10;y2=1-(19*exp(-t2/2)/20;plot(t1,y1,o,t2,y2,k);lege nd(od e,dsolve);年级、专业 2012级数学与应用数学姓名 张旭 学号 12160011063 名单序号16实验一常微分方程实验时间2014年3月19日

12、使用设备、软件PC, MATLAB注：实验报告的最后一部分是实验小结与收获(2) 厂家问：要做多少时间广告，可使市场购买率达到 80% ？4 由解析解可列出方程1 - (19*exp(-t/2)/20=0.8 ，所以t - -2ln -19解得 t=3.11636. (肿瘤生长)肿瘤大小V生长的速率与 V的a次方成正比，其中a为形状参数，0虫乞1; 而其比例系数K随时间减小，减小速率又与当时的 K值成正比，比例系数为环境参数 b。 设某肿瘤参数a=1, b=0.1, K的初始值为2，V的初始值为1。问(1)此肿瘤生长不会超过多大？k =k = -0.1k由已知列出方程组丿，代入具体数值，得丿，

13、v=kvav = kvk(0) =2,v(0) -1，(k,v都是关于时间t的函数)解析法：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,L in eWidth,2); xlabel(t 轴);ylabel(v 轴);因肿瘤不断长大，故t趋于无穷时，该肿瘤达到最大，此时极限为exp(20)=4.8517*10A8，故此肿瘤生长不会超过4.8517*10人89 / 92

14、014春数学实验实验一常微分方程实验时间2014年3月19日使用设备、软件PC, MATLAB注：实验报告的最后一部分是实验小结与收获(2)过多长时间肿瘤大小翻一倍？令 exp(20)*exp(-20*exp(-t/10)=2 ，解得 t=-10*ln(1-1/20*ln2)=0.3527.(3) 何时肿瘤生长速率由递增转为递减？由已求得的结果可得v 与t的关系为v =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); plot(t1,v1,Li neWidth,2);

15、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若参数a=2/3 呢？k=-0.1k1、建立方程组2， k(0) = 2, v(0) =1V = kv3解析法：k,v=dsolve(Dk=-0.1*k,Dv=k*vA(2/3), k(0)=2,v(0)=1,t)k = 2*exp(-t/10)2*exp(-t/10)2*exp(-t/10)v = -(20*exp(-t/10) - 23)A

16、3/27年级、专业 2012级数学与应用数学姓名 张旭 学号 12160011063 名单序号16实验时间2014年3月19 日使用设备、软件PC, MATLAB实验一常微分方程注:实验报告的最后一部分是实验小结与收获(37/2 + (3人(1/2)*3许)/2 - 20*exp(-t/10)A3/27-(20*exp(-t/10) + (3A(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).A3/27;plo

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

18、abel(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.(生态系统的振荡现象)第一次世界大战中，因为战争很少捕鱼，按理战后应能捕到更11 / 92014春数学实验实验一常微分方程注:实验报告的最后一部分是实验小结与收获多的鱼才是。可是大战后，在地中海却捕不到鲨鱼，因而渔民大惑不解。令X1为鱼饵的数量，X2为鲨鱼的数量，t为时间。微分方程为dt式中a1, a 2, b 1, b 2都是正常数。第一式鱼饵X1的增长速度大体上与 X1成正比，即按a1X1比率增加，而被鲨鱼吃掉的部分按 b1X1X2的比率减少；第二式中鲨鱼的增长速度由 于生存竞争的自然死亡或互相咬食按a2X2的比率减少，但又根据鱼饵的量的变化按b2X1 X2 的比率增加。对 a 1=3, b 1=2, a 2=2.5, b2=1, x 1（0）=x