谁帮我看看这个分岔图对不对啊

无水1324

这些都是随便的，没有定量的要求，图能够直观显示就可以了

sssssxxxxx921

时域图是不是也是去掉瞬态解部分？　绘制的时候是不是所取计算周期应该更长一些　

有点不敢确定　呵呵

无水1324

没有必要，相图都稳定了，时域的肯定也稳定了的

sssssxxxxx921

在用数值算法（如龙格库塔法）计算得出：时间tc=[tc,t_span(8000:n)];
位移 yc=[yc,q(8000:n,2)];
我用以下ＭＡＴＬＡＢ语句怎么显示不出来幅值谱图啊　　帮忙看看有什么错误
这些天麻烦你了　　太感谢你了
Y = fft(yc,1024);
Pyy = Y.* conj(Y) / 1024;
f = 1000*(0:1023)/1024;
plot(f(1:512),Pyy(1:512))

无水1324

位移 yc=[yc,q(8000:n,2)];你用这个做什么？

唯位移直接用,yc=q(8000:n,1)];

sssssxxxxx921

是有点费事，呵呵，
不过那和上边的结果没什么关系呀
做不出来幅值谱图　　这是为什么啊

无水1324

你把完整的程序全部发上来我帮你改吧，你这个问题发了很多相关的帖子

sssssxxxxx921

好的　谢谢了　　

sssssxxxxx921

其中有一些参数是做庞加莱图时用到的　没来得及删
这是主程序：（要求２００００转时的时域信号和幅值谱图）
function BallBrg_NonL_Forum
% 求解外圈固定球轴承的变柔度(VC-Varying Compliance)振动（基于赵凌燕的论文）
% 程序有一些不合理、甚至错误的地方，可以用更好的代码代替，由于时间关系没有修改，
% 如有人感兴趣可以把修改的程序发布出来。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 作者：toes
% 版本：论坛发布版
% 相关程序：BallBrg_NonL_Sub_Forum
% 调试环境：Matlab7.0   WinXP SP2
% 参考文献：
% 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc

%% 参数设置
% 用了全局变量来传递一些变量，不推荐，但是懒得改了，好心人优化一下。
global w d D Nb gama kn M C F

% 为了方便绘制分岔图而设置的参数
n_One_T = 100;% 每个周期的采样点数
n_T = 100;% 采样时间占几个周期

% 61903/P5(17*30*7) 球轴承参数
d=0.0173;% 内滚道直径
D=0.0265;% 外滚道直径
Nb=9;     % 滚子数

n_n = 0;
w_limit1=100;% 最低转速（rpm）
w_limit2=20000;% 最高转速（rpm）
w_step = 100;% 转速变化步长(rpm)

q_initial(1:4,1) = 1e-11;% 初始值
gama = 0.00002;% 间隙(m)
F = 6;% 径向力(N)
kn = 7.055e9*0.001^1.5;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.5)。赵的论文给出。
M=0.6*[1 0;0 1];% 质量矩阵
C=200*[1 0; 0 1];% 阻尼矩阵

%% 响应计算
xc=[];
wc=[];
yc=[];xcc=[];ycc=[];tc=[];

    w = 20000*pi/30;% 转化为rad/s单位
   
    wi = w;% 内圈角速度
    wo = 0;% 外圈角速度
   
    w_cage = ( wi*d/2+wo*D/2 )/2/((D+d)/4);% 保持架
    w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
    T_vc = 1/w_vc;% vc周期

    dt=T_vc/n_One_T;% 取点时间步长，dt随转速变化。
    time=n_T*T_vc;% 总的时间
   
    n = round(time/dt);% 离散点数
    t_span(1:n) = linspace(0,time,n);% 时间数组
   
    [t,q]= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
    % 至于用什么ode函数求解合适需要比较验证

  
tc=[tc,t_span(8000:n)];
xc=[xc,q(8000:n,1)];
yc=[yc,q(8000:n,2)];
xcc=[xcc,q(8000:n,3)];
ycc=[ycc,q(8000:n,4)];

Y = fft(yc,1024);
Pyy = Y.* conj(Y) / 1024;
f = 1000*(0:1023)/1024;
plot(f(1:512),Pyy(1:512))

这是子程序：
function dq = BallBrg_NonL_Sub_Forum(t,q)
% BallBrg_NonL调用的微分方程子程序
% 求解外圈固定球轴承的变柔度(VC)振动
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 作者：toes
% 版本：论坛发布版
% 相关程序：BallBrg_NonL_Forum
% 参考文献：
% 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global w d D Nb gama kn M C F

wi = w;
wo = 0;
w_cage=( wi*d+wo*D )/4/((D+d)/4);% 保持架转速（rad/s）

fq=zeros(2,1);% 轴承力初值

diff_1_3 = q(1,1);% 水平方向位移
diff_2_4 = q(2,1);% 垂直方向位移

% 求轴承的非线性反力
for No_ball=1:Nb
    sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
    Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
        + diff_2_4*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
    % 判断哪几个滚动体受到接触力
    if Clearance(No_ball)<=0;
       Clearance(No_ball) = 0;
    end
    fs = abs( (1000*Clearance(No_ball))^1.5 );

    fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball));
    fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
end

F_m1d1_cos = 0;% 不平衡力在水平方向的投影。本例不考虑。
F_m1d1_sin = 0;% 不平衡力在垂直方向的投影。本例不考虑。

Fq(1,1)= - fq(1,1) + F_m1d1_cos;% 水平方向外力
Fq(2,1)= - fq(2,1) + F_m1d1_sin + F;% 垂直方向外力

K = [0 0; 0 0];% 刚性转子，轴段为刚性。

% 动力学微分方程
dq(3:4,1)=inv(M)*(Fq-K*q(1:2,1)-C*q(3:4,1));% x和y方向加速度
dq(1:2,1)=q(3:4,1);

无水1324

我都看蒙呢
佩服你的编程能力，再仔细看一下

sssxxx921

呵呵　　不好意思　见笑了　
这程序主体是别人的　　我只是在他上边改了改　

sssssxxxxx921

上午用我哥的号回复你的　　怎么下午看不到我的贴子了

无水1324

tc=[tc,t_span(8000:n)];
xc=[xc,q(8000:n,1)];
yc=[yc,q(8000:n,2)];
xcc=[xcc,q(8000:n,3)];
ycc=[ycc,q(8000:n,4)];
这个有什么意义？

tc=t_span(8000:n)
xc=q(8000:n,1)
yc=q(8000:n,2)
xcc=q(8000:n,3)
ycc=q(8000:n,4)

sssssxxxxx921

那是我画分岔图和庞加莱图时用的　　改的时候没有去掉

无水1324

这个看了一下，没有什么大的问题，只是有一些语句，你真的要改一下，可以简单明了一点。