2013春数学实验基础 实验报告(1)常微分方程.doc

年级、专业 数学与应用数学 姓名 王锴丰 学号 150142011035 名单序号 实验时间 2013年 3月 18日 使用设备、软件 PC, MATLAB 注： 实验报告的最后一部分是实验小结与收获 实验一 常微分方程1. 分别用Euler法和ode45解下列常微分方程并与解析解比较： (1) 先编写Euler法的M文件function x,y = euler(odefun,tspan,y0,h)x = tspan(1):h:tspan(2);y(1) = y0;for i = 1:length(x)-1 y(i+1) = y(i)+h*feval(odefun,x(i),y(i);endx = x;y = y;1、ode45解 fun = inline(x+y,x,y); ode45(fun,0,3,1)2、 Euler法: t,y=euler(fun,0,3,1,0.01); plot(t,y,r)3、 解析解： s=dsolve(Dy=y+x,y(0)=1,x) s = 2*exp(x) - x - 1 x=0:0.05:3; s=2*exp(x) - x - 1; plot(x,s,:)可以看到三者几乎重合。（2）现将其改成方程组形式：y=y(1);y=y(2);则有 1、ode45解 fun=inline(y(2);0.01*y(2)2-2*y(1)+sin(t),t,y); ode45(fun,0,5,0,1)结果如图二2、解析解： s=dsolve(D2y-0.01*Dy2+2*y-sin(t),y(0)=0,Dy(0)=1,t)Warning: Explicit solution could not be found. In dsolve at 101s = empty sym 解析解在无法解出2. 求一通过原点的曲线，它在处的切线斜率等于若上限增为1.58，1.60会发生什么？ 解：相当于用解微分方程 输入如下指令： fun=inline(2*x+y2,x,y); ode45(fun,0,1.57,0) ode45(fun,0,1.58,0) ode45(fun,0,1.60,0)比较三个图像x上限为1.57时y的上限160左右当x上限为1.58时y的上限达到18乘以10的13次方数量级当x上限为1.60时y的上限还是18乘以10的13次方数量级说明1.57到1.58附近函数开始上升得特别快再增加,解会爆炸3. 求解刚性方程组：解：ode45法： fun=inline(-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5,x,y); tic; x,y=ode45(fun,0,50,1,-1); toc;Elapsed time is 162.804781 seconds.解析解： S=dsolve(Df=-1000.25*f+999.75*g+0.5,Dg=999.75*f-1000.25*g+0.5,f(0)=1,g(0)=-1); S.f,S.g ans = 1/exp(2000*t) - 1/exp(t/2) + 1 ans =1 - 1/exp(2000*t) - 1/exp(t/2)4. (温度过程)夏天把开有空调的室内一支读数为20的温度计放到户外，10分钟后读25.2, 再过10分钟后读数28.32。建立一个较合理的模型来推算户外温度。解：由温度与温差的关系可列微分方程 T=k(c-T),T(0)=20 输入：dsolve(DT=k*(c-T),T(0)=20,t) 得ans = c - (c - 20)/exp(k*t)利用T(10)=25.2, T(20)=28.32拟合(或者解非线性方程) fun=inline(c(1)+exp(-c(2)*t)*(-c(1)+20),c,t) ;lsqcurvefit(fun,30 1,10 20,25.2 28.32)ans = 33.0000 0.0511即解得户外温度c=33,比例系数k=0.05. 5. (广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。(1) 建立该问题的数学模型，并求其数值解与模拟结果作以比较(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？解：(1)由题意可知微分方程 x/x=0.5*(1-x),x(0)=0.05 解析解： dsolve(Dx/x=0.5*(1-x),x(0)=0.05) ans = 1/(exp(log(19) - t/2) + 1) t=0:0.01:10; x=1./(exp(log(19) - t/2) + 1); plot(t,x) hold on数值解： fun=inline(0.5*(1-x)*x,t,x); t,x=ode45(fun,0 10,0.05); plot(t,x,o) %两种解比较如右图(2)、 id=min(find(x0.8); t(id) ans =8.6731即有做8.6731，销售率达到80%6. (肿瘤生长) 肿瘤大小V生长的速率与V的a次方成正比，其中a为形状参数，0a1;而其比例系数K随时间减小，减小速率又与当时的K值成正比，比例系数为环境参数b。设某肿瘤参数a=1, b=0.1, K的初始值为2，V的初始值为1。问(1)此肿瘤生长不会超过多大？(2)过多长时间肿瘤大小翻一倍？(3)何时肿瘤生长速率由递增转为递减？(4)若参数a=2/3呢？解：由题意可知微分方程为：V(t)=K(t)*V(t)a,K(t)=-b*K(t) 令y(1)=V(t),y(2)=K(t); （1）可直接求解析解，再由极限观察可知肿瘤不超过exp(20)的大小（2）数值解微分方程组： fun=inline(y(2)*y(1)1;-0.1*y(2),t,y) ; x,y=ode45(fun,0,1,1,2); plot(x,y) id=min(find(y2.0); x(id) ans = 0.3750 %此时间后肿瘤变为两倍（3）由微分方程可知，肿瘤变化的速率为V(t)则速率变化率为V(t)=K(t)*V(t)+K(t)*V(t)=-b*K(t)*V(t)+K(t)2*V(t)=-b*y(2)*y(1)+y(2)2*y(1) s=y(:,1).*(-0.1).*y(:,2)+y(:,2).2.*y(:,1); it=min(find(s x(it)ans = 30.5892 %即到30单位时间速率开始减小（4）若a=2/3,用同样的方法：得到解析解得到微分方程组的解： S=dsolve(Df=g*f(2/3),Dg=-0.1*g,f(0)=1,g(0)=2); S.f ans = -(20/exp(t/10) - 23)3/27 (37/2 - 20/exp(t/10) + (3*3(1/2)*i)/2)3/27 -(20/exp(t/10) - 37/2 + (3*3(1/2)*i)/2)3/27再求极限： a=limit(S.f,t,inf) a = 12167/27 %此为实部极限值为450.6296 (37/2 + (3*3(1/2)*i)/2)3/27 -(- 37/2 + (3*3(1/2)*i)/2)3/27故肿瘤最大不能超过451单位体积,又解得在t=0.4时达到两倍 t= 9.5718时速率开始减小选做题：1.(生态系统的振荡现象)第一次世界大战中，因为战争很少捕鱼，按理战后应能捕到更多的鱼才是。可是大战后，在地中海却捕不到鲨鱼，因而渔民大惑不解。令x1为鱼饵的数量，x2为鲨鱼的数量，t为时间。微分方程为 (5.20)式中a1, a2, b1, b2都是正常数。第一式鱼饵x1的增长速度大体上与x1成正比，即按a1x1比率增加, 而被鲨鱼吃掉的部分按b1x1x2的比率减少；第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a2x2的比率减少，但又根据鱼饵的量的变化按b2x1x2的比率增加。对a1=3, b1=2, a2=2.5, b2=1, x1(0)=x2(0)=1求解。画出解曲线图和相轨线图，可以观察到鱼饵和鲨鱼数量的周期振荡现象。解：将a1=3, b1=2, a2=2.5, b2=1, x1(0)=x2(0)=1代入求解微分方程组。 fun=inline(x(1)*(3-2*x(2),-x(2)*(2.5-x(1),t,x); x,t=ode45(fun,0,10,1,1); plot(t,x(:,2),:) hold on plot(t,x(:,1),r) plot(x(:,1),x(:,2),r)解曲线图如右上：可以看到鲨鱼和鱼饵的数量呈大致相同的周期变化，周期T约为2.5相轨线图如右下：可观察到鲨鱼数量和鱼饵数量基本维持一闭合曲线的方程呈周期变化。实验小结： 本次实验通过对有关的常微分方程的