请高手解答~~解非线性油膜力模型

lhy

m文件如下：
function dx=ymwd(t,x)
global w prou;
prou=0.1;
%w=7600;
m=2;    %kg
C=0.1; %mm
n=0.018;%Pa.s
R=25;%mm
L=10;%mm
g=9.8;%N/kg^2   %omg=sqrt(m*qs/(m*c));
qs=n*w*R*L*(R/C)^2*(L/2/R)^2;
M=(w^2*m*C)/qs;
G1=g/C/w^2;
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1).^2-x(3).^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1).^2-x(3).^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a)).^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3).^2);
fx=sqrt((x(1)-2*x(4)).^2+(x(3)+2*x(2)).^2)/(1-x(1).^2-x(3).^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4)).^2+(x(3)+2*x(2)).^2)/(1-x(1).^2-x(3).^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
dx=[x(2);fx/M+prou*cos(t);x(4);fy/M+prou*sin(t)-G1];
用龙格库塔方法解出来的时域波形为一条斜线，如下图：哪里出错了？请高手指点~~

octopussheng

这是刚性的一种体现，需要调整你的系统参数再计算了。

lhy

回复 2 # octopussheng 的帖子

请问具体到哪些系统参数的改变呢？