计算物理解 常微分方程编程.doc

例1：function jyjsT,Y=ode45(yjs,0:10,2,1);plot(T,Y(:,1),T,Y(:,2)function ydot=yjs(t,y)ydot=y(2); 4 ;例2：研究受空气阻尼的抛体运动function znxpglobal m d e %用全局变量传递参数m=1;b=0,0.2,0.2;p=0,0,1; %传递给方程的参数px=4.6;4.5;4.5; % 说明文字的x坐标py=3.5;1.8;0.4; % 说明文字的y坐标strdd1=无阻尼; %图形说明文字strdd2=阻力正比于v; strdd3=阻力正比于v2;figurefor i=1:3 %三次循环计算三种条件下的运动d=b(i); e=p(i); %选定传递参数的值t,y=ode45(znxpfun,0:0.01:10,0,3,0,5);Hi=max(y(:,3) %轨迹最高点Ti=t(find(y(:,3)=Hi) %到达最高点时间vx0i=y(find(y(:,3)=Hi),2) %到最高点速度subplot(2,1,1) %分区画图axis(0 6 -70 2); % 确定坐标范围hold on %保留图形xlabel(x); ylabel(y); %x,y轴的标注plot(y(:,1),y(:,3); %粒子的空间轨迹subplot(2,1,2) axis(0 10 0 4) hold on xlabel(t);ylabel(dx/dt)text(px(i),py(i),strddi); %加注解文字plot(t,y(:,2) %速度随时间变化图end%-function ydot=znxpfun(t,y) %要解的常微分方程global m d e ydot=y(2); -d/m*y(2)*(y(2).2+y(4).2)(e/2); y(4); -9.8-d/m*y(4)*(y(2).2+y(4).2)(e/2); 例3: 混沌理论的洛仑兹方程 function lor%建立函数文件axis(10 50 -50 50 -50 50) %设定坐标轴范围view(3) %设定视角hold on %保持设置title(Lorenz Attractor) %设定标题options =odeset(OutputFcn,odephas3) ; %设定选项,绘出解向量中前三个元素序列构成的三维相空间图t,y=ode23(lorfun,0,20,0,0,eps,options); %解方程sol=ode23(lorfun,0,20,0,0,eps) %以结构数组输出解ss=deval(sol,10) %查看t=10时的解function ydot=lorfun(t,y)ydot=-8/3*y(1)+y(2)*y(3); -10*y(2)+10*y(3); -y(2)*y(1)+35*y(2)-y(3); 例4：小球弹跳ballode 例5：求解常微分方程F=(t,y)y(2);-y(1);%F=(t,y)0,1;-1,0*y; %使用矩阵形式来编写匿名函数%ode45(F,0,10,2,0)sol=ode45(F,0,10,0,1)ss=deval(sol,5.15)options=odeset(outputfcn,odephas2);%改变原来输出的设置ode45(F,0,10,0,1,options)例5：火焰问题（刚性问题）delta=0.0001;F=inline(y2-y3,t,y);opts=odeset(RelTol,1.e-4);subplot(1,2,1)ode45(F,0,2/delta,delta,opts);subplot(1,2,2)ode23s(F,0,2/delta,delta,opts);%用符号工具箱求解y=dsolve(Dy=y2-y3,y(0)=0.01);y=simplify(y);pretty(y) %以书写习惯显示符号表达式的函数ezplot(y,0,200)例6：事件功能eventsfunction shijianopts=odeset(events,g);y0=1;0;t,y,tfinal=ode45(f,0,Inf,y0,opts);tfinalplot(t,y(:,1),-,0,tfinal,1 0,o)function ydot=f(t,y) %模型函数文件ydot=y(2); -1+y(2)2;function gstop,isterminal,direction=g(t,y) %检测事件发生的函数文件 gstop=y(1); isterminal=1; direction=; 例7：弹簧摆运动function thbglobal L m k gtheta0=pi/10; %初始角度,可由读者设不同的值m=1;k=80;g=9.8; L0=1;L=L0+m*g/k; %L0为弹簧原来长度，L为弹簧静止时长度t,u1=ode45(thbfun,0:0.005:15,L0 0 theta0 0);y1,x1=pol2cart(u1(:,3),u1(:,1);y1 = -y1;%将极坐标换为直角坐标 figureymax=max(abs(y1);axis(-1.2 1.2 -1.2*ymax 0.2);%设置坐标范围axis offtitle(弹簧摆,fontsize,14)hold on; R =0.055 ; %设置弹簧半径yy = -L0 : 0.01 : 0;xx = R*sin(yy./L0*30*pi);%用正弦曲线表示弹簧a,r = cart2pol(xx,yy); %用坐标变换来画初始位置的弹簧a = a + theta0;xx,yy = pol2cart(a,r);%弹簧的数据 line(-1 1,0 0,color,r,linewidth,2)%画横杆，球与弹簧ball = line(xx(1),yy(1),color,r,marker,.,markersize,70,erasemode,xor);ball2 = line(xx(1),yy(1),color,0.5 0.51 0.6,linestyle,-,linewidth,1.3,erasemode,none);spring = line(xx,yy,color,r,linewidth,2,erasemode, xor);pause(0.5) for i = 1 : length(t) %画出每一步的弹簧位置 yy = -u1(i,1) : 0.01 : 0; xx = R*sin(yy./u1(i,1)*30*pi); a,r = cart2pol(xx,yy); a = a + u1(i,3); xx,yy = pol2cart(a,r); set(ball,XData,x1(i),YData,y1(i); set(ball2,XData,x1(i),YData,y1(i); set(spring,XData,xx,YData,yy); drawnow;end %子函数文件function F = thbfun(t,u)global L m k gF = u(2); u(1)*u(4)2 + g*cos(u(3) - k/m*(u(1) - L + m*g/k); u(4); -2*u(2)*u(4)/u(1) - g*sin(u(3)/u(1);例8：圆锥陀螺运动%程序1按照常规写法解圆锥陀螺运动function tlydR=1;h=2;p=1;g=9.8;%陀螺参数a=2/(4*h2/R2+1); b=5*p*h/(4*h2+R2);phi=0;thi=pi/6;psi=0; %初始条件%options=odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5 1e-4 1e-4 1e-5);t,u=ode45(tlydfun,0:0.02:10,thi 0 phi 0 psi 75); figure;axis(-1 1 -1 1 0 1);hold on ;text(-0.5,2,1,陀螺运动,fontsize,20,color,r);grid on; x0,y0,z0=cylinder(0:0.05:0.5,60);%陀螺侧面的数据网络Q=linspace(0,2*pi,60);%以下四句是陀螺上表面面的数据网络cx0=0.5*cos(Q);cy0=0.5*sin(Q);cz0=ones(1,60);x,y,z=zbbh(x0,y0,z0,thi,phi,psi);%陀螺初始位置侧面数据cx,cy,cz=zbbh(cx0,cy0,cz0,thi,phi,psi);%陀螺初始位置上表面数据xm0=0;ym0=0;zm0=1.2; %中轴线数据xm,ym,zm=zbbh(xm0,ym0,zm0,thi,phi,psi); %中轴线h3=plot3(0,xm,0,ym,0,zm,r,linewidth,5);h1=surf(x,y,z);%画初始位置陀螺侧面h2=fill3(cx,cy,cz,m);%画初始位置陀螺上表面h4=plot3(cx(1),cy(1),cz(1),b.,markersize,20);%画初始位置陀螺轴pause(0.5); for i=1 :length(u); x,y,z=zbbh(x0,y0,z0,u(i,1),u(i,3),u(i,5); cx,cy,cz=zbbh(cx0,cy0,cz0,u(i,1),u(i,3),u(i,5); cn=xm,ym,zm; xm,ym,zm=zbbh(xm0,ym0,zm0,u(i,1),u(i,3),u(i,5); set(h1,xdata,x,ydata,y,zdata,z); set(h2,xdata,cx,ydata,cy,zdata,cz); set(h3,xdata,0,xm,ydata,0,ym,zdata,0,zm); set(h4,xdata,cx(1),ydata,cy(1),zdata,cz(1); drawnow; end text(-2,0,0,结束,fontsize,30,color,r);%计算加速度u(:,7)=u(:,4).2.*sin(u(:,1).*cos(u(:,1)*(1-a)-a*u(:,4).*u(:,6).*sin(u(:,1)+b.*g.*sin(u(:,1);u(:,8)=u(:,2).*u(:,4).*cot(u(:,1).*(a-2)+a.*u(:,2).*u(:,6)./sin(u(:,1);u(:,9)=u(:,2).*u(:,4).*(sin(u(:,1)-cot(u(:,1).*cos(u(:,1).*(a-2)-a.*u(:,2).*u(:,6).*cot(u(:,1);text_s=%章动角-t,进动角-t,自转角-t,%章动角速度-t,进动角速度-t,自转角速度-t,. %章动角加速度-t,进动角加速度-t,自转角加速度-t,;figure(Units,normalized,Position,0.05 0.1 0.9 0.8,. Tag,curve,NumberTitle,off,Name,%欧拉角特征曲线);for i=1:3 for j=1:3 axes(Position,0.05+0.32*(j-1) 0.04+0.32*(3-i) 0.27 0.24,. XGrid,on,YGrid,on); uicontrol(Units,normalized,Position,0.11+0.32*(j-1) 0.28+0.32*(3-i) 0.15 0.04,. Style,text,BackgroundColor,0.8 0.8 0.8,. String,text_s(i-1)*3+j,FontSize,12); plot(t,u(:,(i-1)*3+j); endend%-坐标变换函数-function x,y,z=zbbh(x0,y0,z0,thi,phi,psi)x=x0*(cos(psi)*cos(phi)-sin(psi)*cos(thi)*sin(phi). +y0*(-sin(psi)*cos(phi)-cos(psi)*cos(thi)*sin(phi). +z0*sin(thi)*sin(phi);y=x0*(cos(psi)*sin(phi)+sin(psi)*cos(thi)*cos(phi). +y0*(-sin(psi)*sin(phi)+cos(psi)*cos(thi)*cos(phi). -z0*sin(thi)*cos(phi);z=x0*sin(psi)*sin(thi)+y0*cos(psi)*sin(thi)+z0*cos(thi);%-需求解的微分方程组-function y = tlydfun(t,u)R=1;h=2;p=1;g=9.8;a=2/(4*h2/R2+1); b=5*p*h/(4*h2+R2);y=u(2); (1-a)*u(4)2*sin(u(1)*cos(u(1)-a*u(4)*u(6)*sin(u(1)+b*g*sin(u(1); u(4); (a-2)*u(2)*u(4)*cot(u(1)+a*u(2)*u(6)/sin(u(1); u(6); u(2)*u(4)*(sin(u(1)-(a-2)*cot(u(1)*cos(u(1)-a*u(2)*u(6)*cot(u(1); %程序2书上圆锥陀螺运动程序function tlydclcclear allR=1;h=2;p=1;g=9.8;%陀螺参数,R为圆锥底面半径, h为高, p为陀螺材质的密度, g为重力加速度a=2/(4*h2/R2+1); b=5*p*h/(4*h2+R2);phi=0;thi=pi/6;psi=0; %初始条件t,u=ode45(tlydfun, 0:0.02:10,thi phi psi 0 0 75); axis(-1 1 -1 1 0 1);grid on;title(%拉格朗日-泊松情况下的圆锥陀螺运动)hold on ;%text(-0.5,2,1,陀螺运动,fontsize,20,color,r);x0,y0,z0=cylinder(0:0.05:0.5,60);%陀螺侧面的数据网络Q=linspace(0,2*pi,60);%以下四句是陀螺上表面面的数据网络cx0=0.5*cos(Q);cy0=0.5*sin(Q);cz0=ones(1,60);x,y,z=zbbh(x0,y0,z0,thi,phi,psi);%陀螺初始位置侧面数据cx,cy,cz=zbbh(cx0,cy0,cz0,thi,phi,psi);%陀螺初始位置上表面数据xm0=0;ym0=0;zm0=1.2; %中轴线数据xm,ym,zm=zbbh(xm0,ym0,zm0,thi,phi,psi); %中轴线h3=plot3(0,xm,0,ym,0,zm,r,linewidth,5);h1=surf(x,y,z);%画初始位置陀螺侧面h2=fill3(cx,cy,cz,m);%画初始位置陀螺上表面h4=plot3(cx(1),cy(1),cz(1),b.,markersize,20); %陀螺面起点的位置，每个时刻不同pause(0.5);%在继续执行之前,暂停0.5秒 for i=1 :length(u); x,y,z=zbbh(x0,y0,z0,u(i,1),u(i,2),u(i,3); cx,cy,cz=zbbh(cx0,cy0,cz0,u(i,1),u(i,2),u(i,3); %求陀螺新位置的数据 cn=xm,ym,zm; xm,ym,zm=zbbh(xm0,ym0,zm0,u(i,1),u(i,2),u(i,3); %中轴线 set(h1,xdata,x,ydata,y,zdata,z); %画新位置的陀螺 set(h2,xdata,cx,ydata,cy,zdata,cz); set(h3,xdata,0,xm,ydata,0,ym,zdata,0,zm); set(h4,xdata,cx(1),ydata,cy(1),zdata,cz(1); plot3(cn(1),xm,cn(2),ym,cn(3),zm); %留下中轴轨迹 drawnow; end %计算加速度u(:,7)=u(:,5).2.*sin(u(:,1).*cos(u(:,1)*(1-a)-a*u(:,5).*u(:,6).*sin(u(:,1)+b.*g.*sin(u(:,1);u(:,8)=u(:,4).*u(:,5).*cot(u(:,1).*(a-2)+a.*u(:,4).*u(:,6)./sin(u(:,1);u(:,9)=u(:,4).*u(:,5).*(sin(u(:,1)-cot(u(:,1).*cos(u(:,1).*(a-2)-a.*u(:,4).*u(:,6).*cot(u(:,1);text_s=%章动角-t,进动角-t,自转角-t,%章动角速度-t,进动角速度-t,自转角速度-t,. %章动角加速度-t,进动角加速度-t,自转角加速度-t,;figure(Units,normalized,Position,0.05 0.1 0.9 0.8,. Tag,curve,NumberTitle,