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

1. 分别用Euler法和ode45解下列常微分方程并与解析解比较： (1) 编写Euler法的matlab函数，命名为euler.mfunction t,y=euler(odefun,tspan,y0,h)t=tspan(1):h:tspan(2);y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h*feval(odefun,t(i),y(i);endt=t;y=y;下面比较三者的差别：% ode45odefun=inline(x+y,x,y);x1,y1=ode45(odefun,0,3,1);plot(x1,y1,ko);pausehold on;% Eulerx2,y2=euler(odefun,0,3,1,0.05);plot(x2,y2,r+);pausehold on;% 解析解y0=dsolve(Dy=t+y,y(0)=1);ezplot(y0,0,3);pausehold off;legend(ode45,euler法,解析解);Euler法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。 (2) 先写M文件ex1_2fun.mfunction f=ex1_2fun(t,y)f(1)=y(2);f(2)=0.01*y(2).2-2*y(1)+sin(t);f=f(:);% ode45t1,y1=ode45(ex1_2fun,0,5,0;1);plot(t1,y1(:,1),ko);% 解析解s=dsolve(D2y-0.01*(Dy)2+2*y=sin(t),y(0)=0,Dy(0)=1,t)s = empty sym %由此可知该微分方程无解析解2. 求一通过原点的曲线，它在处的切线斜率等于若上限增为1.58，1.60会发生什么？ odefun=inline(2*x+y2,x,y);subplot(1,4,1);x1,y1=ode45(odefun,0,1.57,0);plot(x1,y1,r*);title(上限1.57);subplot(1,4,2);x2,y2=ode45(odefun,0,1.58,0);plot(x2,y2,bo);title(上限1.58);subplot(1,4,3);x3,y3=ode45(odefun,0,1.6,0);plot(x3,y3,k);title(上限1.60);subplot(1,4,4);plot(x1,y1,r*);hold on;plot(x2,y2,bo);hold on;plot(x3,y3,k);hold off;legend(上限1.57,上限1.58,上限1.60);结论：随着x上界的增加，解趋于无穷大。3. 求解刚性方程组：先写M函数ex3fun.mfunction f=ex3fun(t,y)f(1)=-1000.25*y(1)+999.75*y(2)+0.5;f(2)=999.75*y(1)-1000.25*y(2)+0.5;f=f(:);%作图t,y=ode15s(ex3fun,0,50,1,-1);plot(t,y,*);4. (温度过程)夏天把开有空调的室内一支读数为20的温度计放到户外，10分钟后读25.2, 再过10分钟后读数28.32。建立一个较合理的模型来推算户外温度。设：t时刻温度计的读数为T，户外温度为c，T的增速与室内外温差（c-T）成正比，由此建立微分方程，其中k为比例系数%首先，计算解析解 y=dsolve(DT=k*(c-T),T(0)=20,t)y =c - (c - 20)/exp(k*t)%又已知，用非线性最小二乘拟合该函数，调用lsqcurvefit命令：fun=inline(c(1)-(c(1)-20)./exp(c(2)*t),c,t);lsqcurvefit(fun,30 1,10 20,25.2 28.32)ans = 33.0000 0.0511即户外温度c=33，比例系数k=0.05115. (广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。(1) 建立该问题的数学模型，并求其数值解与模拟结果作以比较；设：市场占有率为，市场占有率的增长速度为，则相对增长率为，由此建立微分方程为%首先，计算解析解 s=dsolve(Dx=(0.5*(1-x)*x,x(0)=0.05,t)s=1/(exp(log(19) - t/2) + 1)%再调用ode45计算数值解，并作图比较解析解与数值解的区别：odefun=inline(0.5*(1-x)*x,t,x);t,x=ode45(odefun,0,20,0.05);plot(t,x,r*);hold on;ezplot(s,0 20);hold off;(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？t_min=min(find(x0.8);t(t_min)ans = 8.8543结果：大约8.8543个单位的时间后，可使市场购买率达到80%。6. (肿瘤生长) 肿瘤大小V生长的速率与V的a次方成正比，其中a为形状参数，;而其比例系数K随时间减小，减小速率又与当时的K值成正比，比例系数为环境参数b。设某肿瘤参数a=1, b=0.1, K的初始值为2，V的初始值为1。问(1)此肿瘤生长不会超过多大？由已知条件建立微分方程：%先编写上述函数的ex6_1fun.m文件function f=ex6_1fun(t,y)f(1)=y(2).*y(1);f(2)=-0.1*y(2);f=f(:);%画出其图像，并求最大的肿瘤大小Vt1,y1=ode45(ex6_1fun,0,100,1,2);plot(t1,y1(:,1),r*);max(y1(:,1)ans = 4.8567e+008%故肿瘤生长不会超过(2)过多长时间肿瘤大小翻一倍？%从图像可以看出，大约在（0，1）内，肿瘤大小翻一倍，以此求解t2,y2=ode45(ex6_1fun,0 1,1;2);t_min=min(find(y2(:,1)2);t2(t_min)ans =0.3750答：大约0.4个单位时间后肿瘤大小翻一倍。(3)何时肿瘤生长速率由递增转为递减？dv=y1(:,2).*y1(:,1);Vn,tn=max(dv);t1(tn)ans = 29.5466答：大约30个单位时间后肿瘤生长速率由递增转为递减。(4)若参数a=2/3呢？%重新编写微分方程ex6_4fun，并依次计算上述三个问题function f=ex6_4fun(t,y)f(1)=y(2)*y(1).(2/3);f(2)=-0.1*y(2);f=f(:);%画图像，肿瘤生长不会超过450.7959%求多长时间肿瘤大小翻一倍，大约为0.4ans = 0.4000%何时肿瘤生长速率由递增转为递减，大约为10ans = 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求解。画出解曲线图和相轨线图，可以观察到鱼饵和鲨鱼数量的周期振荡现象。%首次编写上述微分方程的ex_7fun.m函数function f=ex_7fun(t,x)a1=3;b1=2;a2=2.5;b2=1;f(1)=x(1)*(a1-b1*x(2);f(2)=-x(2)*(a2-b2*x(1);f=f(:);%画出解曲线图和相轨线图t,x=ode45(ex_7fun,0,4,1;1);subplot(1,2,1);plot(t,x(:,1),r,t,x(:,2),k:);legend(x1,x2);title(解曲线);subplot(1,2,2);plot(x(:,1),x(:,2);title(相轨线);2.编写四阶Runge-Kutta法程