大虾帮助，如何用matlab画以下方程的极限环

byandby

dz=[z(2);-c1/m1*z(2)-k1/m1*z(1);z(4);-c2/m2*z(4)-k2/m2*z(3)]+[0;Ff/m1;0;-Ff/m2];
m1=m2=1;k1=k2=1,c1=c1=0.01,
mus=0.6;
N=10;alpha=0.012;v0=1;
Pn=mus*N;
Ff0=N*(mus-alpha*v0);
Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
if abs(vr)==0
    Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>0
    Ff=(mus-alpha*vr)*N*sign(vr)-Ff0;
end
if -vr>0
    Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end

无水1324

vr 和x是什么？是否固定的还是随响应变化的？

另外数值求解已经分析可以做出来吗？

byandby

vr是变化的，是两物体之间的相对速度，x是初始值，x=[0 0 0 0;0 0.01 0 0.01];我自己编了一个如下：
function dz=nianhuamoban(t,z,c1m1,k1m1,Ffm1,c2m2,k2m2,Ffm2)
dz=[z(2);-c1m1*z(2)-k1m1*z(1);z(4);-c2m2*z(4)-k2m2*z(3)]+[0;Ffm1;0;-Ffm2];

clear;
clc;
ep=0.005;
m1=1;m2=1;k1=1;k2=1;c1=0.01;c2=0.01;mus=0.6;
N=15;alpha=0.012;v0=1;
Pn=mus*N;
Ff0=N*(mus-alpha*v0);
c1m1=c1/m1;k1m1=k1/m1;c2m2=c2/m2;k2m2=k2/m2;
a=1;                    
b=5;
x=[0 0 0 0;0 0.01 0 1];     %质量刚度相等时，两个速度相等时或位移相等时收敛
Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
vr=v0+x(2,4)-x(2,2);      
if abs(vr)<=ep
    Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>ep
    Ff=(mus-alpha*vr)*N*sign(vr)-Ff0;
end
if -vr>ep
    Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end
Fs1(1)=Fs;
h=1;
y(1,:)=x(2,:);
Ffm1=Ff/m1;
Ffm2=Ff/m2;
t0(1)=a*0.1;
while(1)
       h=h+1;
       tspan=(a:1:b)*0.1;
       [t,x]=ode45(@nianhuamoban,tspan,x(2,:),[],c1m1,k1m1,Ffm1,c2m2,k2m2,Ffm2);
       Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
       Fs1(h)=Fs;
       vr=v0+x(2,4)-x(2,2);
       y(h,:)=x(2,:);
       if abs(vr)<ep
          Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
       end
       if vr>=ep
          Ff=(mus-alpha*vr)*N*sign(vr)-Ff0;
       end
       if -vr>=ep
          Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
       end  
%        g(h)=y(h,2)-y(h-1,2);
%        if abs(g(h))<=0.001
%            y(h,2)=y(h-1,2);
%        end
%        G(h)=y(h,4)-y(h-1,4);
%        if abs(G(h))<=0.001
%            y(h,4)=y(h-1,4);
%        end
       Ffm1=Ff/m1;
       Ffm2=Ff/m2;
       t0(h)=t(2);
      
       a=a+1;
       b=b+1;
       if a==10000
           break;
       end
end
t00=t0';
% y=y/w1*u0;
figure(1)
subplot(2,1,1);
plot(t00,y(:,1));
subplot(2,1,2);
plot(t00,y(:,2));
figure(2)
plot(y(:,1),y(:,2));
figure(3)
subplot(2,1,1);
plot(t00,y(:,3));
subplot(2,1,2);
plot(t00,y(:,4));
figure(4)
plot(y(:,3),y(:,4));
但是出来结果跟文献里的不相符啊，麻烦无水前辈指教

[ 本帖最后由 byandby 于 2009-12-12 16:47 编辑 ]

无水1324

Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
vr=v0+x(2,4)-x(2,2);      
if abs(vr)<=ep
    Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>ep
    Ff=(mus-alpha*vr)*N*sign(vr)-Ff0;
end
if -vr>ep
    Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end

应该直接放在nianhuamoban函数里面

byandby

谢谢无水前辈耐心解答，您的意思是那时候已经得到X的值了吗？按照您这么说的话就不用循环了，是吧？

[ 本帖最后由 byandby 于 2009-12-12 20:11 编辑 ]

无水1324

需要，但是你要知道每一个x值与Ff有关系啊 ，所以这个是要放在里面的。他是方程的一个参数

byandby

非常感谢前辈的解答，可我循环中不是已经有Ffm1=Ff/m1; Ffm2=Ff/m2了，Ffm1和Ffm2也代入下一次循环的呀

无水1324

但是必须在方程的m函数里面哈

byandby

是啊，无水前辈，我定义了两个函数，分别是nianhuamoban和solvenianhuamoban，我自己觉得循环没问题啊，可是结果中并没有出现极限环的粘滑现象啊，无水前辈能否帮忙分析分析啊

byandby

我定义的两个函数如下，要达到的结果如下

无水1324

难道你真的不能把循环放在m程序里面去？

byandby

但是按照无水无水前辈的意思，怎么定义x这个数组变量？x首先是初值，然后又是每一次的迭代结果。怎么赋值呢？

byandby

按照无水前辈的改法：
function dz=nianhuamoban(t,z)
ep=0.005;
m1=1;m2=1;k1=1;k2=1;c1=0.01;c2=0.01;mus=0.6;
N=15;alpha=0.012;v0=1;
Pn=mus*N;
Ff0=N*(mus-alpha*v0);
c1m1=c1/m1;k1m1=k1/m1;c2m2=c2/m2;k2m2=k2/m2;
Fs=k1*z(1)+c1*z(2)-k2*z(3)-c2*z(4);
vr=v0+z(4)-z(2);      
if abs(vr)<=ep
    Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>ep
    Ff=(mus-alpha*vr)*N*sign(vr)-Ff0;
end
if -vr>ep
    Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end
Ffm1=Ff/m1;
Ffm2=Ff/m2;
dz=[z(2);-c1m1*z(2)-k1m1*z(1);z(4);-c2m2*z(4)-k2m2*z(3)]+[0;Ffm1;0;-Ffm2];

主程序：
clear;
clc;
m1=1;m2=1;k1=1;k2=1;c1=0.01;c2=0.01;mus=0.6;
a=1;                    
b=5;
x=[0 0 0 0;0 0.01 0 1];     %质量刚度相等时，两个速度相等时或位移相等时收敛
Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
Fs1(1)=Fs;
h=1;
y(1,:)=x(2,:);
t0(1)=a*0.1;
while(1)
       h=h+1;
       tspan=(a:1:b)*0.1;
       [t,x]=ode45(@nianhuamoban,tspan,x(2,:));
       Fs1(h)=Fs;
       y(h,:)=x(2,:);      
       t0(h)=t(2);      
       a=a+1;
       b=b+1;
       if a==1000
           break;
       end
end
t00=t0';
figure(1)
subplot(2,1,1);
plot(t00,y(:,1));
subplot(2,1,2);
plot(t00,y(:,2));
figure(2)
plot(y(:,1),y(:,2));
figure(3)
subplot(2,1,1);
plot(t00,y(:,3));
subplot(2,1,2);
plot(t00,y(:,4));
figure(4)
plot(y(:,3),y(:,4));
出来的效果见下图，跟没改之前的效果一样啊，还是出不了文献中的效果，请无水前辈赐教。还有，极限环就是根据相图画的吧

无水1324

能不能把这篇文章拿来看看

byandby

前辈客气了，当然可以啊