求助分岔图程序

chuandong418

调了很久了，还是调不出来，发上来，希望能得到大家的帮忙，谢谢了！
主程序：function xdot=111(t,x,m,k,kc,c,f,e,s,jianxi,w0s,w0,g)
xdot=zeros(4,1);
xdot(1)=x(3);
xdot(2)=x(4);
xdot(3)=-(((sqrt(x(1)^2+x(2)^2)-1)/sqrt(x(1)^2+x(2)^2))*kc*(x(1)-f*x(2))...
        *(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2)...
        /(m*w0s)-(c/(m*w0))*x(3)-(k/(m*w0s))*x(1)...
        +(e*s^2*cos(s*t))/jianxi;
           %[-2*kc*x(3)-x(1)+(-(e-jianxi)*k/e)*(x(1)-f*x(2))/(m*w0s)+e*wn^2*cos(w*t)];  
xdot(4)=-(((sqrt(x(1)^2+x(2)^2)-1)/sqrt(x(1)^2+x(2)^2))*kc*(f*x(1)+x(2))...
        *(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2)...
        /(m*w0s)-(c/(m*w0))*x(4)-(k/(m*w0s))*x(2)...
        +(e*s^2*sin(s*t))/jianxi-g/(w0s*jianxi);

运行程序：
clear all
tic
m=4       %      20;  %4.0;                 %kg
k=0.25*10^6       %    6*10^6;
kc=6.0*10^7       %4*10^7                %Nm-1
c=1200;                %Nsm-1
f=0.2    %0.06;
e=0.00006    %0.05;          %m
jianxi=0.00015          %1*10^(-3);   %0.169;    %设置间隙值
w0s=k/m;              %w0平方
w0=sqrt(w0s);          %求w0
kcc=c/(2*m*w0);        %阻尼比
g=9.8;
X0=[0.001;0.001;0;0]; %[0,0.1,2.1,0.2];
rate=[0.1:0.1:7];
offset=1;
i=1;
for i=1:length(rate)
G=rate(i);
tspan=[0 140*pi/rate(i)];
T=[0:pi/(4*rate(i)):140*pi/rate(i)];
w=rate(i)*w0;
[T,Y]=ode45(@111,T,X0',[],m,k,kc,c,f,e,rate(i),jianxi,w0s,w0,g);
t=T';
cutindexes=find(T>100*pi/rate(i));
T=T(cutindexes);
Ytmp1=Y(cutindexes,1);
Ytmp2=Y(cutindexes,2);
Ytmp1(:,2)=Ytmp2(:,1);
Y=Ytmp1;
t=t(cutindexes);
tmp=length(Y);
m(offset:(tmp+offset-1),1)=g(i);
m(offset:(tmp+offset-1),2)=Y(:,2);
offset=tmp+offset;
end
figure
plot(m(:,1),m(:,2),'k.','markersize',4)
title('Bifurcation Diagramm')
xlabel('s')
ylabel('Y')
toc;

老出现那几个错误，望过来人指点，小女子在此谢谢了！

gghhjj

很简单的一个错误

matlab中文件名或者函数名第一个符号不能是数字，只能是字母

这是第一个问题

gghhjj

另外(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2)一段应该是有问题的

如果sqrt中取为负值，则sign得到的值将是复数

[ 本帖最后由 mjhzhjg 于 2007-8-12 16:38 编辑 ]

无水1324

原帖由 chuandong418 于 2007-5-23 22:47 发表

                               
登录/注册后可看大图

调了很久了，还是调不出来，发上来，希望能得到大家的帮忙，谢谢了！
主程序：function xdot=111(t,x,m,k,kc,c,f,e,s,jianxi,w0s,w0,g)
xdot=zeros(4,1);
xdot(1)=x(3);
xdot(2)=x(4);
xdot(3)=-(((sqrt(x ...

还有就是编辑主函数的时候
function xdot=***(t,x,m,k,kc,c,f,e,s,jianxi,w0s,w0,g)
应该改为：
function xdot=***(t,x,flag,m,k,kc,c,f,e,s,jianxi,w0s,w0,g)

***为含有字母的函数名。

[ 本帖最后由 无水1324 于 2007-5-24 10:30 编辑 ]

chuandong418

55555,按两位学长提到的去修改了，但是还是有错误！痛苦呢！
??? Error: File: TESTXdot.m Line: 1 Column: 28
Incomplete or misformed expression or statement.

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

咕噜噜

xdot(3)=-(((sqrt(x(1)^2+x(2)^2)-1)/sqrt(x(1)^2+x(2)^2))*kc*(x(1)-f*x(2))...
        *(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2)...
        /(m*w0s)-(c/(m*w0))*x(3)-(k/(m*w0s))*x(1)...
        +(e*s^2*cos(s*t))/jianxi;
我这里显示这里也有错误，但是总也改不好:@L

咕噜噜

clear all
tic
w0s=k/m;              %w0平方
w0=sqrt(w0s);          %求w0
kcc=c/(2*m*w0);        %阻尼比
g=9.8;
X0=[0.001;0.001;0;0]; %[0,0.1,2.1,0.2];
rate=[0.1:0.1:7];
offset=1;
for i=1:length(rate)
G=rate(i);
tspan=[0 140*pi/rate(i)];
T=[0:pi/(4*rate(i)):140*pi/rate(i)];
w=rate(i)*w0;
[T,Y]=ode45(@fun,T,X0',[],m,k,kc,c,f,e,G,jianxi,w0s,w0,g);
t=T';
cutindexes=find(T>100*pi/rate(i));
T=T(cutindexes);
Ytmp1=Y(cutindexes,1);
Ytmp2=Y(cutindexes,2);
Ytmp1(:,2)=Ytmp2(:,1);
Y=Ytmp1;
t=t(cutindexes);
tmp=length(Y);
m(offset:(tmp+offset-1),1)=g;
m(offset:(tmp+offset-1),2)=Y(:,2);
offset=tmp+offset;
end
figure
plot(m(:,1),m(:,2),'k.','markersize',4)
title('Bifurcation Diagramm')
xlabel('s')
ylabel('Y')
toc;


function xdot=fun(t,x,m,k,kc,c,f,e,s,jianxi,w0s,w0,g)
m=4;       %      20;  %4.0;                 %kg
k=0.25*10^6;       %    6*10^6;
kc=6.0*10^7;       %4*10^7                %Nm-1
c=1200;                %Nsm-1
f=0.2;    %0.06;
e=0.00006;    %0.05;          %m
jianxi=0.00015;  
xdot=zeros(4,1);
xdot(1)=x(3);
xdot(2)=x(4);
xdot(3)=-(((sqrt(x(1)^2+x(2)^2)-1)/sqrt(x(1)^2+x(2)^2))*kc*(x(1)-f*x(2)) ...
        *(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2) ...
        /(m*w0s)-(c/(m*w0))*x(3)-(k/(m*w0s))*x(1) ...
        +(e*s^2*cos(s*t))/jianxi;
           %[-2*kc*x(3)-x(1)+(-(e-jianxi)*k/e)*(x(1)-f*x(2))/(m*w0s)+e*wn^2*cos(w*t)];  
xdot(4)=-(((sqrt(x(1)^2+x(2)^2)-1)/sqrt(x(1)^2+x(2)^2))*kc*(f*x(1)+x(2)) ...
        *(sign(abs(sqrt(x(1)^2+x(2)^2)-1))+sign(sqrt(x(1)^2+x(2)^2)-1))/2) ...
        /(m*w0s)-(c/(m*w0))*x(4)-(k/(m*w0s))*x(2) ...
        +(e*s^2*sin(s*t))/jianxi-g/(w0s*jianxi);
这样调出来时没问题的，但是画出来的图形就比较怪了，因为m(offset:(tmp+offset-1),1)=g;也就是矩阵m的第一列都等于一个常数，但是你原来的g(i)的形式肯定是不对的，因为g这里是一个常量，或者你是不是写错了？

gghhjj

最好把方程给出来，用word或者图片的形式
方便大家对照你的程序

chuandong418

上面程序的方程如下：

gghhjj

看了一下你的方程，建议仔细点，感觉方程描述明显不对
F和Fx、Fy好像没关系，另外ks如何确定？

doori99

你的问题解决了吗?
能否给我发个程序
我也刚开始做转子碰摩
谢谢
邮箱
doori99@126.com

[ 本帖最后由 无水1324 于 2007-8-7 22:15 编辑 ]