实验一常微分方程

《实验一常微分方程》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、实用文档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);y(1)=y0;fori=1:length(x)-1y(k+1)=y(i)+h*feval(odefun,x(i),y(i));endx=x';y=y';endode45:odefun=inline('x+y','x','y');xspan=[0,3];实用文档y0=1

2、;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',x3,y3,'k*');xlabel('x轴');ylabel('y轴');legend('euler','ode45','dsolve');ode45求得的结果与用解析法求得的结果更接近，故ode45的精度较高，Euler法求得的结果精度较低。(2)令则原方程等价于方程组：，，不能解析，只能用数值法求解。Euler：实用文档

3、function[t,y]=euler2(odefun1,odefun2,tspan,y0,h)t=tspan(1):h:tspan(2);y(1,1)=y0(1);y(2,1)=y0(2);fori=1:length(t)-1k1=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;endt=t';y=y';endode45：odefun1=inline('0*t1+0*y1+y2');odefun2=inline('-2*

4、y1+0.01*y2^2+sin(t1)');[t1,y1]=euler2(odefun1,odefun2,[0,5],[0,1],0.1);[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：functiondy=eu(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-2*y(1)+0.01*y(2)^2+sin(t);od

5、e45求得的结果精度较高，euler法求得的结果在准确值上下波动。2.一通过原点的曲线，它在处的切线斜率等于若上限增为1.58，1.60会发生什么？等价于求解，且的初值问题。解析法：y=dsolve('Dy=2*x+y^2','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+y^2');实用文档subplot(1,3,1);o

6、de45(odefun,[0,1.57],0);title('01.5之后的斜率增长速度很快，若上限增为1.58，1.60，则相应的y将会出现更大的增长。3.求解刚性方程组：functiondy=fun(t,y)dy=zeros(2,1);dy(1)=-1000.25*y(1)+999.75*y(2)+0

7、.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℃。建立一个较合理的模型来推算户外温度。由题意可知由于随着时间的增加，温度增长越来越慢，户外温度与温度计