声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6833|回复: 18

[稳定性与分岔] 大虾帮助,如何用matlab画以下方程的极限环

[复制链接]
发表于 2009-12-12 10:38 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
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
回复
分享到:

使用道具 举报

发表于 2009-12-12 15:47 | 显示全部楼层
vr 和x是什么?是否固定的还是随响应变化的?

另外数值求解已经分析可以做出来吗?
 楼主| 发表于 2009-12-12 16:29 | 显示全部楼层
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 编辑 ]

想要的效果

想要的效果
发表于 2009-12-12 16:52 | 显示全部楼层
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函数里面
 楼主| 发表于 2009-12-12 20:08 | 显示全部楼层
谢谢无水前辈耐心解答,您的意思是那时候已经得到X的值了吗?按照您这么说的话就不用循环了,是吧?

[ 本帖最后由 byandby 于 2009-12-12 20:11 编辑 ]
发表于 2009-12-12 22:15 | 显示全部楼层

回复 5楼 byandby 的帖子

需要,但是你要知道每一个x值与Ff有关系啊 ,所以这个是要放在里面的。他是方程的一个参数
 楼主| 发表于 2009-12-12 22:28 | 显示全部楼层
非常感谢前辈的解答,可我循环中不是已经有Ffm1=Ff/m1; Ffm2=Ff/m2了,Ffm1和Ffm2也代入下一次循环的呀
发表于 2009-12-16 16:51 | 显示全部楼层
但是必须在方程的m函数里面哈
 楼主| 发表于 2009-12-17 22:05 | 显示全部楼层

回复 8楼 无水1324 的帖子

是啊,无水前辈,我定义了两个函数,分别是nianhuamoban和solvenianhuamoban,我自己觉得循环没问题啊,可是结果中并没有出现极限环的粘滑现象啊,无水前辈能否帮忙分析分析啊
 楼主| 发表于 2009-12-17 22:13 | 显示全部楼层
我定义的两个函数如下,要达到的结果如下

函数1

函数1

函数2

函数2

函数2

函数2

函数2

函数2

想达到的结果1

想达到的结果1

想达到的结果1

想达到的结果1
发表于 2009-12-18 08:08 | 显示全部楼层
难道你真的不能把循环放在m程序里面去?
 楼主| 发表于 2009-12-19 14:15 | 显示全部楼层

回复 11楼 无水1324 的帖子

但是按照无水无水前辈的意思,怎么定义x这个数组变量?x首先是初值,然后又是每一次的迭代结果。怎么赋值呢?
 楼主| 发表于 2009-12-19 14:46 | 显示全部楼层
按照无水前辈的改法:
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));
出来的效果见下图,跟没改之前的效果一样啊,还是出不了文献中的效果,请无水前辈赐教。还有,极限环就是根据相图画的吧

相图1

相图1

相图2

相图2
发表于 2009-12-24 09:46 | 显示全部楼层

回复 13楼 byandby 的帖子

能不能把这篇文章拿来看看
 楼主| 发表于 2009-12-24 19:14 | 显示全部楼层
前辈客气了,当然可以啊

1

1

2

2

3

3

4

4

5

5

6

6

7

7

8

8

9

9

10

10

11

11

12

12
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-17 11:53 , Processed in 0.070447 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表