这是我的变系数动力学方程： 按照书中308页下边的那种龙格库塔的方法，结合您的帖子上面的程序，不做变换化为一阶方程，而直接用龙哥-库塔对其求解。我的编程思路大致如下： clc clear %%%%%%%%%%%%%%%%%%%%%%% Ixy=0.81*10^-8; Ixz=0.56*10^-2; Iyx=0.81*10^-8; Iyy=0.66432; Iyz=1.166*10^-7; Izx=0.56*10^-2; Izy=1.166*10^-7; Izz=0.64; 密度，质量 转动惯量等已知是的赋值 %%%%%%%%%%%%%%%%%%%%%%%%%%%% tt=0:0.01:1; %%%%%%这里也试过0.001,0.0001等 for ij=1:length(tt) t=tt(ij) dt=0.01;%%%%%%%%%步长也试过0.1,0.0001等 。。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。 M= K= F= C= 。。。。。。。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。。。。。%%%%%句号中间是求M K F C （都是关于t的函数） Inv_M=inv(M); if ij==1 q=zeros(112,1);%%%%初始位移 q1d=zeros(112,1); %%%%初始速度 q2d=Inv_M*(QA-K*q-C*q1d); %%%%初始加速度 end if ij>1 qq=q; qq1d=q1d; qq2d=q2d; Kd1=Inv_M*(-K*qq-C*qq1d+F);% Kd2=Inv_M*(-K*(qq+0.5*dt*qq1d)-C*(qq1d+0.5*dt*Kd1)+F); Kd3=Inv_M*(-K*(qq+0.5*dt*qq1d+0.25*dt^2*Kd1)-C*(qq1d+0.5*dt*Kd2)+F); Kd4=Inv_M*(-K*(qq+dt*qq1d+0.5*dt^2*Kd2)-C*(qq1d+dt*Kd3)+F); update_disp=qq+dt*qq1d+(dt^2)/6*(Kd1+Kd2+Kd3);% update_velo=qq1d+dt/6*(Kd1+2*Kd2+2*Kd3+Kd4); update_acc =Kd1; q=update_disp; q1d=update_velo; q2d=update_acc; q107=q(107); %%%x位移 q108=q(108); %%%y位移 q109=q(109); %%%y位移 qqq1(ij)=q107; qqq2(ij)=q108; qqq3(ij)=q109; end end figure(1) plot(tt,qqq1,'b')%%%x位移 title('λòÆ') xlabel('ê±¼ät/s') ylabel('λòÆ/m') figure(2) plot(tt,qqq2,'m')%%%y位移 title('Ëù¶è') xlabel('ê±¼ät/s') ylabel('λòÆ/m') figure(3) plot(tt,qqq3,'m')%%%z位移 title('Ëù¶è') xlabel('ê±¼ät/s') ylabel('λòÆ/m')
| 
发表于 2015-10-14 09:49
显身卡
发表于 2015-10-14 11:46
楼主
