常微分方程及其求解实例资料课件.ppt

《常微分方程及其求解实例资料课件.ppt》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、第4章常微分方程及其求解常微分方程（组）数值解1初值问题一阶常微分方程一阶常微分方程高阶常微分方程(1)Euler法function[tout,yout]=euler(ypfun,t0,tfinal,y0,tol,trace)pow=1/3;ifnargin<5,tol=1e-3;endifnargin<6,trace=0;endt=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;tout=zeros(chunk,1);yout=zeros(chunk,length(y));k

2、=1;tout(k)=t;yout(k,:)=y.';iftraceclc,t,h,yend用差商近似代替导数while(tt)ift+h>tfinal,h=tfinal-t;endf=feval(ypfun,t,y);f=f(:);delta=norm(h*f,'inf');tau=tol*max(norm(y,'inf'),1.0);ifdelta<=taut=t+h;y=y+h*f;k=k+1;ifk>length(tout)tout=[tout;zeros(chunk,1)];yout=[y

3、out;zeros(chunk,length(y))];endtout(k)=t;yout(k,:)=y.';endiftracehome,t,h,yendifdelta~=0.0h=min(hmax,0.9*h*(tau/delta)^pow);endendif(t

4、('ypfun',0,1,1)Euler法用差商近似代替导数(2)Runge-Kutta法梯形公式展开二阶Runge-Kutta法辛普森公式展开由Euler法估计三阶Runge-Kutta法二阶、三阶Runge-Kutta函数ode23[t,y]=ode23(@fun,tspan,y0)[t,y]=ode23(@fun,tspan,y0,options,p1,p2)四阶、五阶Runge-Kutta函数ode45[t,y]=ode45(@fun,tspan,y0)[t,y]=ode45(@fun,tspan,y0,options

5、,p1,p2)求微分方程初始条件[DyDt.m]functionydot=DyDt(t,y)mu=2;ydot=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];spin=[0,30];y0=[1;0];[tt,yy]=ode45(@DyDt,tspan,y0);plot(tt,yy(:,1))xlabel('t'),title('x(t)')plot(yy(:,1),yy(:,2))xlabel('位移'),ylabel('速度')围绕地球旋转的卫星轨道问题的形成轨道上运动的卫星，在Newton第二定律，和万有引

6、力定律作用下假定卫星以初速度处进入轨道。在构成一阶微分方程组[dYdt.m]functionYd=DYdt(t,Y)%t一定是标量形式的自变量%Y必须是列向量globalGME%在函数中定义全局变量传递参数xy=Y(1:2);Vxy=Y(3:4);%前两个元素是"位移量"，后两个是"速度量"。r=sqrt(sum(xy.^2));DYdt1=Vxy(1);DYdt2=Vxy(2);DYdt3=-G*ME*xy(1)/r^3;DYdt4=-G*ME*xy(2)/r^3;Yd=[DYdt1;DYdt2;DYdt3;DYdt4];%

7、Yd必须是与Y同维的列向量。globalGME%在主程序中定义全局变量传递参数G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=[t0,tf];%指定解算微分方程的时间区间Y0=[x0;0;0;vy0];%给定初值向量[t,YY]=ode45('DYDt',tspan,Y0);X=YY(:,1);%输出Y的第一列是位移数据x(t)Y=YY(:,2);%输出Y的第二列是位移数据y(t)plot(X,Y,'b','Linewidth',2);hol

8、donaxis('image')%保证x、y轴等长刻度，且坐标框恰包容图形[XE,YE,ZE]=sphere(10);%产生单位球面数据RE=0.64e7;%地球半径XE=RE*XE;YE=RE*YE;ZE=0*ZE;%坐标纸上的地球平面数据mesh(XE,YE,ZE),ho