该楼层疑似违规已被系统折叠隐藏此楼查看此楼
如下是我写的计算步骤，稍后解释：(不想看可以跳过~)
%main
%input some parameters
m1=input('please input the mass of your leg: ');
m2=input('please input the mass of your upper body: ');
h1=0.5*input('please input the length of your leg: ');
h2=0.5*input('pleae input the length of your upper body:');
sita=input('please input the angle of dip: ');
ga=input('please input an angle that reflects your sensitivity on reflection: (should in (0,pi))');
r=input('please input the radius of the aparatus:');
m0=input('please input the mass of the aparatus:');
g=9.8;
%transfer matrix
A=[cos(sita) 0 sin(sita);0 1 0;-sin(sita) 0 cos(sita)];
%angle between your leg and your body
beta=(2*sita/pi)*(alpha-ga);
%centroid vector in the 'double index' frame
cm11=[r*cos(alpha) r*sin(alpha) 0];
cm22=[(r+h2*sin(beta))*cos(alpha) (r+h2*sin(beta))*sin(alpha) 0];
%centroid position in the normal frame
cm1=(A*cm11')';
cm2=(A*cm22')';
%moment of inertia
J=m1*r^2+0.5*m0*r^2+m2*((4/3)*(h2^2)*(sin(beta))^2+r^2+2*r*h2*sin(beta));
%unit vector of the rotation axis in both frames
n0=[0 0 1];
n=(A*n0')';
%moment to the axis due to gravity
component=dot(n0,n);
gra1=m1.*g.*n0-component.*n;
gra2=m2.*g.*n0-component.*n;
moment1=cross(cm1,gra1);
moment2=cross(cm2,gra2);
%dynamic differential equation
J*(D2alpha)=dot(n,moment1)+dot(n,moment2);
大意就是一个力学模型，现在需要解一个二阶常微分方程但是很复杂，大概是有(D2y=f(y)的形式，f是一个由一些向量点乘、叉乘blabla得出来的以y为自变量的函数，但主要是三角函数与多项式、分式函数迭代，或者说，这个模型中角度是时间的函数，而惯量张量与力矩都与角度有关= =)。。。。我没有办法化简，求问这种情况怎么求数值解啊！！！
