请大家看一下我的ode45计算是否有问题

chunshui2003

%%%%%方程建立%%%%%%



function uu=fangcheng_x1(t,u);

%%基本参数
m1=170;
m2=90;
e=0.0001;
r=0.0006;
k2=1154744;
omega=19*2*pi;
omega0=15*2*pi;
phi0=0.6/180*pi;          %%%转换为pi
theta0=0.4/180*pi;
xi=0.14;



%%合成参数
phi=omega*t+phi0;         
theta=omega*t+theta0;

M=m1+m2;
A1=m1*e*omega^2*cos(phi);
A2=m2*r*omega.^2*cos(theta);
A3=k2*r*cos(theta);
B1=m1*e*omega.^2*sin(phi);
B2=m2*r*omega.^2*sin(theta);
B3=k2*r*sin(theta);

%%%  未进行参数简化之前的方程
% uu=zeros(4,1);
% uu(1)=u(2);
% uu(2)=(m1*e*omega.^2*cos(omega*t+phi0)+m2*r*omega.^2*cos(omega*t+theta0)- k2*r*cos(omega*t+theta0))./(m1+m2)-2*xi*omega0*u(2)-omega0.^2*u(1);
% uu(3)=u(4);
% uu(4)=(m1*e*omega.^2*sin(omega*t+phi0)+m2*r*omega.^2*sin(omega*t+theta0)- k2*r*sin(omega*t+theta0))./(m1+m2)-2*xi*omega0*u(4)-omega0.^2*u(3);
%其中u1=x=x1,u2=dx/dt,u3=y=y1,u4=dy/dt



%%%% 进行参数简化后的方程
uu=zeros(4,1);
uu(1)=u(2);
uu(2)=(A1+A2-A3)./M-2*xi*omega0*u(2)-omega0.^2*u(1);
uu(3)=u(4);
uu(4)=(B1+B2-B3)./M-2*xi*omega0*u(4)-omega0.^2*u(3);


%%%%% 求解%%%%%

clc
clear
period=1/19;     %%%  旋转频率为19Hz，则计算周期为1/19
tspan=[0:period/100:100*period];
y0=[1 1 1 1];    %% 计算初值
[t,u]=ode45(@fangcheng_x1,tspan,y0);




我觉得我建立的方程和求解没有问题，但是结果和文献上的相差甚远，也直接导致我自己研究内容的后续工作进展不顺利，因为总是怀疑算法出了问题。请无水，oct，咕噜噜等各位朋友们帮忙看一下，谢谢了。

zdeming

有空时我来算算

NASA

是不是文献的单位有问题
圆频率的单位怎么是HZ？

kangarooli

你的肯定有问题啊，时间历程是收敛的，轨迹怎么可能是一条线呢，应该是稳定焦点形式啊！
看你是想画稳定极限环的形式，而得到的是稳定焦点，如果方程没错，那就是你的参数可能有问题，这个问题是很烦人啊！我之前也是，仿真做不出来，直接没法继续下面的工作

zdeming

LZ 不知如何解决的。
该论文的题目是？想查一下。觉得该论文是否有点问题啊？
w0=sqrt((k1+k2)/(m1+m2)是否正确，我知道对于齿轮系统的固有频率 w0
w0=sqrt(k/me)，me为等效质量,me=m1*m2/(m1+m2)，另外：
kesi也有点问题。。