1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1)30,1)0(,<<=+='x y y x y 编写Euler 法的matlab 函数，命名为euler.m function [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));end t=t';y=y';下面比较三者的差别：% ode45odefun=inline('x+y','x','y');[x1,y1]=ode45(odefun,[0,3],1);plot(x1,y1,'ko');pause hold on ;% Euler·¨[x2,y2]=euler(odefun,[0,3],1,0.05);plot(x2,y2,'r+');pause hold on ;% 解析解y0=dsolve('Dy=t+y','y(0)=1');ezplot(y0,[0,3]);pause hold off ;legend('ode45','euler 法','解析解');Euler 法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。

(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<<先写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(:);% ode45[t1,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会(,)x y 22,0 1.57.x y x +<<x 发生什么？odefun=inline('2*x+y^2','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');上上1.5713上上1.5814上上1.6014结论：随着x 上界的增加，解趋于无穷大。

3. 求解刚性方程组：500,1)0(,5.025.100075.999,1)0(,5.075.99925.100022121211<<⎩⎨⎧-=+-='=++-='x y y y y y y y y 先写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 为比例系数⎩⎨⎧=-='20)0()(T T c k T %首先，计算解析解>> y=dsolve('DT=k*(c-T)','T(0)=20','t')y =c - (c - 20)/exp(k*t)%又已知，用非线性最小二乘拟合该函数，调用lsqcurvefit 命令：32.28)20(,2.25)10(==T T 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) 建立该问题的数学模型，并求其数值解与模拟结果作以比较；设：市场占有率为，市场占有率的增长速度为，则相对增长率为，由此建立微分方程为x x 'xx '⎪⎩⎪⎨⎧=-⨯='05.0)0()1(5.0x x x x %首先，计算解析解>> 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 ;t1/(exp(log(19) - t/2) + 1)(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？t_min=min(find(x>0.8));t(t_min)ans = 8.8543结果：大约8.8543个单位的时间后，可使市场购买率达到80%。

6. (肿瘤生长) 肿瘤大小V 生长的速率与V 的a 次方成正比，其中a 为形状参数，;而其比例系数10≤≤a K 随时间减小，减小速率又与当时的K 值成正比，比例系数为环境参数b 。

设某肿瘤参数a=1, b=0.1, K 的初始值为2，V 的初始值为1。

问(1)此肿瘤生长不会超过多大？由已知条件建立微分方程：⎪⎪⎩⎪⎪⎨⎧==⨯-='≤≤⨯='2)0(1)0()()(10,)()()(K V t K b t K a t V t K t V a %先编写上述函数的ex6_1fun.m 文件function f=ex6_1fun(t,y)f(1)=y(2).*y(1);f(2)=-0.1*y(2);f=f(:);%画出其图像，并求最大的肿瘤大小V[t1,y1]=ode45(@ex6_1fun,[0,100],[1,2]);plot(t1,y1(:,1),'r*');max(y1(:,1))8ans =4.8567e+008%故肿瘤生长不会超过8108567.4⨯(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法程序实验小结与收获：完成第二份实验报告后，我学习了Matlab有关求解常微分方程（组）的指令，了解了什么是Euler法和刚性方程组以及如何求解，并了解了ode45、ode15s、Euler等函数间的异同；某些无法求得解析解的方程可以通过计算数值解了解它们的性质，有助于我们的学习；还学习了通过建立数学模型解决实际问题，将理论应用于实际，锻炼了我们逻辑思维能力和抽象概括能力。