mntlyd.txt

一个关于陀螺转动的简单模拟程序,是用MATLAB编写的.

回答网友xiaoxiaoxinghu的程序[一个关于陀螺转动的程序]:

function tlyd;
% tlyd陀螺运动的微分方程

Dfun=inline(['[y(2);(1-a)*(y(4))^2*sin(y(1))*cos(y(1))',...
    '-a*y(4)*y(6)*sin(y(1))+b*g*sin(y(1));',...
    'y(4);(a-2)*y(2)*y(4)*cot(y(1))+a*y(2)*y(6)*csc(y(1));',...
    'y(6);y(2)*y(4)*[sin(y(1))-(a-2)*cot(y(1))*cos(y(1))]',...
    '-a*y(2)*y(6)*cot(y(1))]'],'t','y','flag','a','b','g');
% [t,Y]=ode45(Dfun,[0,20],Y0,[],a,b,g);
% plot(t,Y)

axis([-1 1 -1 1 0 1]);
hold on;
text(0,2,0.8,'陀螺运动','fontsize',40,'color','r');
grid on;
[x0,y0,z0]=cylinder(0:.05:.5,60); %获得陀螺侧面的三维数据网格(固连坐标x0,y0,z0为常数)
axis equal
Q=linspace(0,2*pi,60);%%以下四句是陀螺底面的数据网格(固连坐标cx0,cy0,cz0为常数)
cx0=0.5*cos(Q);
cy0=0.5*sin(Q);
cz0=ones(1,60);
phi=0;thi=pi/6;psi=0;
[x,y,z]=zbbh(x0,y0,z0,thi,phi,psi);%%陀螺初始位置侧面数据(静坐标x,y,z随陀螺的运动而变化)
[cx,cy,cz]=zbbh(cx0,cy0,cz0,thi,phi,psi);%%陀螺初始位置底面数据(静坐标cx,cy,cz随陀螺的运动而变化)
h1=surf(x,y,z);%%画初始位置陀螺侧面并获取所画图形的图柄
h2=fill3(cx,cy,cz,[0 1 1]);%%画初始位置陀螺的底面并获取所画图形的图柄
pause(0.5);%在继续执行之前,暂停0.5秒
R=1;h=2;p=1;g=9.8;%%陀螺的参数, R为圆锥底面半径, h为高, p为陀螺材质的密度, g为重力加速度
a=2/(4*h^2/R^2+1); b=5*p*h/(4*h^2+R^2);
[t,u]=ode45(Dfun,[0:0.1:25],[thi;0;phi;0;psi;75],[ ],a,b,g);%用函数句柄@tlydfun来调用描述陀螺运动微分方程的函数tlydfun
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));%%求陀螺新位置底面数据
    set(h1,'xdata',x,'ydata',y,'zdata',z);%%画新位置陀螺侧面
    set(h2,'xdata',cx,'ydata',cy,'zdata',cz);%%画新位置陀螺底面
    drawnow;
end;
text(-1.4,0,'结束','fontsize',40,'color','r');

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);