实验一 常微分方程.doc

1. 分别用Euler法和ode45解下列常微分方程并与解析解比较： (1) 解析： y=dsolve(Dy=x+y,y(0)=1,x) y = 2*exp(x) - x - 1 Euler: function [x,y]=euler(odefun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for i=1:length(x)-1 y(k+1)=y(i)+h*feval(odefun,x(i),y(i)); end x=x ; y=y ; end ode45: odefun=inline(x+y,x,y); xspan=[0,3]; y0=1; 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： 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); for i=1:length(t)-1 k1=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; end t=t; y=y; end ode45： odefun1=inline(0*t1+0*y1+y2); odefun2=inline(-2*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： function dy=(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=-2*y(1)+0.01*y(2)^2+sin(t); ode45求得的结果精度较高，euler法求得的结果在准确值上下波动。 一通过原点的曲线，它在处的切线斜率等于若上限增为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);ode45(odefun,[0,1.57],0);title(0x1.57); subplot(1,3,2);ode45(odefun,[0,1.58],0);title(0x1.58); subplot(1,3,3);ode45(odefun,[0,1.60],0);title(0x1.60); 曲线单调递增，且在x1.5之后的斜率增长速度很快，若上限增为1.58，1.60，则相应的y将会出现更大的增长。 3. 求解刚性方程组： function dy=fun(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]=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℃