该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

如下是我写的计算步骤，稍后解释：(不想看可以跳过~)

%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为自变量的函数，但主要是三角函数与多项式、分式函数迭代，或者说，这个模型中角度是时间的函数，而惯量张量与力矩都与角度有关= =)。。。。我没有办法化简，求问这种情况怎么求数值解啊！！！