我这个非线性方程程序怎么求不出来啊 帮忙看看

octopussheng

把你的程序再贴一下吧，我帮你算算试试！

sssssxxxxx921

子程序：
function dq=Rotors_System_Sub_Func(t,q) %???
global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2   %????
r01=0.00002;%m ?????
%r02=0.00005;%m
m1=7.86;%kg     ??
m2=11.93;%kg
W1=m1*9.8;%N       ??
W2=m2*9.8;%N
Kb1=4.6346e10;%N/m^1.5       %??
Kb2=3.4953e9;%N/m^1.1
F=1485;
F1 = 1485/2;
delta_i=1.2443*w2^2*1e-11;
delta_e=2.2183*w1^2*1e-11;
Dm=(117+133.3)/2;
wc=(w2*(1-8/Dm)+(1+8/Dm)*w1)/2;
delta_o=4.912*wc^2*1e-14;
delta_F=8.10*F1^0.925/7.8^0.85*1e-8;
r02=0.00005-(delta_i-delta_o-delta_e-delta_F);
%q                   %???????Q??
Fx11=0;
Fy11=0;   %??????
Fx22=0;
Fy22=0;
for i=1:Nb1
   
    sita(i)=2*pi/Nb1*(i-1)+BN/Nb1*w1*2*pi*t;        %??
    Deformation(i,1)=q(1,1)*cos(sita(i))+q(2,1)*sin(sita(i))-r01;  %??
    if Deformation(i)<=0
        Deformation(i)=0;   %?????????
    end
   
    fx11=Kb1*(Deformation(i))^1.5*cos(sita(i));
    fy11=Kb1*(Deformation(i))^1.5*sin(sita(i));        %?
   
    Fx11=Fx11+fx11;      %??
    Fy11=Fy11+fy11;
   
    %sita2(i)=2*pi/Nb1*(i-1)+BN/Nb1*w2*2*pi*t;
    %Deformation2(i,1)=q(3,1)*cos(sita2(i))+q(4,1)*sin(sita2(i))-r01;
    %if Deformation2(i)<=0
    %    Deformation2(i)=0;
    %end
   
    %fx22=Kb1*(Deformation2(i))^1.5*cos(sita2(i));
   % fy22=Kb1*(Deformation2(i))^1.5*sin(sita2(i));
   
   % Fx22=Fx22+fx22;
   % Fy22=Fy22+fy22;
end
Fx12=Fx11;
Fy12=Fy11;
%Deformation               %?????  ????
%???????
Fx21=0;
Fy21=0;
for j=1:Nb2
    sita(j)=2*pi/Nb2*(j-1)+(Ro2/(Ro2+Ri2)*w1+Ri2/(Ro2+Ri2)*w2)*2*pi*t;
    Deformation1(j,1)=q(3)*cos(sita(j))+q(4)*sin(sita(j))-r02;
    if Deformation1(j)<=0
       Deformation1(j)=0;
    end
   
    fx21=Kb2*(Deformation1(j))^1.1*cos(sita(j));
    fy21=Kb2*(Deformation1(j))^1.1*sin(sita(j));
   
Fx21=Fx21+fx21;
Fy21=Fy21+fy21;
end

%Fx22=0;
%Fy22=0;
for k=1:Nb1
   sita(k)=2*pi/Nb1*(k-1)+BN/Nb1*w2*2*pi*t;
    Deformation2(k,1)=q(3,1)*cos(sita(k))+q(4,1)*sin(sita(k))-r01;
    if Deformation2(k)<=0
       Deformation2(k)=0;
   end
   
    fx22=Kb1*(Deformation2(k))^1.5*cos(sita(k));
    fy22=Kb1*(Deformation2(k))^1.5*sin(sita(k));
   
    Fx22=Fx22+fx22;
    Fy22=Fy22+fy22;
end

P=1470;
Cx11=1500;Cy11=1500;         %???
Cx12=2500;Cy12=2500;
Cx21=2000;Cy21=2000;
Cx22=7000;Cy22=7000;
%dq(1:4,1)=q(5:8,1);
dq(5:8,1)=[-1/m1*((Cx11+Cx12)*q(5,1)+Fx11+Fx12-(W1+Fx21+Cx21*q(7,1)));...
         -1/m1*((Cy11+Cy12)*q(6,1)+Fy11+Fy12-(Fy21+Cy21*q(8,1)));...
          -1/m2*((Cx21+Cx22)*q(7,1)+Fx21+Fx22-(P+W2));...
         -1/m2*((Cy21+Cy22)*q(8,1)+Fy21+Fy22)];
dq(1:4,1)=q(5:8,1);            %????

主程序：
function Rotors_System_Func          %???
clear
clc
global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2         %????
n_one_T=200;           %???????
n_T=100;               %?????
      
Ro1=4.8828e-2;         %??
Ri1=6.8672e-2;
Nb1=10;         %???
  
Ro2=66.515e-3;
Ri2=58.5e-3;            %??????
Nb2=34;        %???
n_n=0;             %????
w2_min=100;       % ??
w2_max=10600;
w2_step=100;         %??????
w1_rpm=11000;      %W1??
q_initial(1:8,1)=1e-11;    %?????
BN=Ri1/(Ri1+Ro1)*Nb1;      %???????
tic
%w2_rpm=100;     
for w2_rpm=w2_min:w2_step:w2_max        
    n_n=n_n+1
    w2=w2_rpm/60;
    w1=w1_rpm/60;               %?????
    r=Ro2/(Ro2+Ri2);
   
    %w_cage_1=w1*r+w2*(1-r);         %?????
    %w_vc_1=w_cage_1*Nb2;
    %T_vc_1=2*pi/w_vc_1;
   
    w_cage=w2*BN;
    w_vc=w_cage;
   
    T_vc=2*pi/w_vc*n_n;        %??
    dt=T_vc/n_one_T;       %????
   
    %n_1=T_vc_1/dt;
   
    time=n_T*T_vc;            %???
    n=round(time/dt);        %????
    t_span(1:n)=linspace(0,time,n);   
    %options=odeset('RelTol',1e-6);
   
    %[t,q]=ode45('Rotors_System_Sub_Func',t_span,q_initial,options);
   [t,q]=ode45('Rotors_System_Sub_Func',t_span,q_initial);
    %subplot(2,2,1);plot(w2_rpm,q(8500:100:10000,5),'k.','markersize',5)
    %hold on
    plot(w2_rpm,q(15000:200:end,8),'k.','markersize',5)
    hold on
    %subplot(2,2,3);plot(w2_rpm,q(8500:100:10000,7),'k.','markersize',1)
    %hold on
    %subplot(2,2,4);plot(w2_rpm,q(8500:100:10000,8),'k.','markersize',1)
    %hold on
end
toc

sssssxxxxx921

上面没有任何注释   
重新发个附件   添加了一点注释

花如月

哇塞，楼修的这么高呀，我也出点力，在精神上支持一下

sssssxxxxx921

谢谢了   现在除了有人陪伴  就很需要精神上的鼓励的了   要不就跳楼了  :lol

octopussheng

我先用你的程序算算，看看结果如何吧！

sssssxxxxx921

计算应该是没问题的就是在选择庞加莱截面时没法选择   三个频率对应的周期  而且是不可通约的
第二个图比第一个图取得周期（准确说计算时间段）要大，得出来图比较好，说明周期第一个图在一个系统周期上取了很多点，与一个周期取一个点（即用庞加莱截面截 解曲线）相矛盾
所以   只有找到系统周期   才能用庞加莱截面上找映射点（即每周期取一个点）

octopussheng

结果出来了，看看吧，好像和你上面的也差不多啊！

sssssxxxxx921

就是一样的

octopussheng

还有没有什么地方值得修改的，你看看，我也帮你算算，呵呵！

sssssxxxxx921

就是在选周期那存在问题    可是多频激励周期按传统的最小公倍数   对于这个不可通约的无理小数行不通    即使固定转速画庞加莱图时，对各个周期做近似  也不能取的好的周期（那样周期太长）   计算点和周期数没法进行

sssssxxxxx921

这是我做的相图      能从这两个图上直接判断周期  拟周期 和混沌吗

sssssxxxxx921

说错了刚才发上去的是轴心轨迹图，
再发几个相图    难道没有一个是稳定解的吗  :@L

咕噜噜

轴心的相图来看1000转是拟周期，4100时看上去也很象拟周期，不知道你取了多少个点？
下面哪几个相图感觉看不出走向，不好说

sssssxxxxx921

再看看这个相图    怎么会这样呢   想不通   第一个图为x1和dx1  第二个图y1和dy1   第三个x2和dx2   第四个图y2和dy2

怎么会出现第二个图那样呢   一个是周期一个是非周期的样子  怎么回事啊