如何利用一下程序作出poincare圖或分叉圖

是滑動軸承在某一參數(shù)下的動力學(xué)特性：

function df=rhs(t,f,data)

%方程中的系統(tǒng)參數(shù)

m=data(1);g=data(2);c=data(3);n=data(4);b=data(5);mf=data(6);M=data(7);D0=data(8);De=data(9);

tau0=data(10);Ke=data(11);K0=data(12);omega=data(13);e=data(14);

%化簡后的動力學(xué)方程

df(1)=f(2);

df(2)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(2)/(M*omega)-2*tau0*mf*f(4)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(1)/(M*omega^2)-(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(3)/(M*omega)+m*e*cos(t)/(M*c);

df(3)=f(4);

df(4)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(4)/(M*omega)+2*tau0*mf*f(2)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(3)/(M*omega^2)+(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(1)/(M*omega)+m*e*sin(t)/(M*c)-m*g/(M*c*omega^2);

df=[df(1);df(2);df(3);df(4)];

運行以下程序

clc;

clear;

format long

m=50;? ?? ? Ke=7.2762e6;? ? De=2e3;

n0=0.079;? ? m0=-0.25;? ? g=9.8;

e=0.0002;? ? R=0.067;? ?? ?l=0.102;

n=2.5;? ?? ? b=0.45;? ???tau0=0.4;? ?c=0.0003;? ?? ?nu=50;

z=0.1;? ?? ? p=2e5;? ?? ???delta=0.0003;? ? T=l/nu;? ???upsilon=1.47e-5;

fid=fopen('poincare.txt','w');

fidd=fopen('time.txt','w');

fiddd=fopen('center.txt','w');

% 轉(zhuǎn)速

omega=1300;

% 初始位置

x0=0;

y0=-m*g/Ke;

X0=x0/c;Y0=y0/c;

% 無量綱轉(zhuǎn)速

nq=length(omega);

for i=1:nq

s(i)=omega(i)/sqrt(Ke/m);

% Child短密封軸承動力系數(shù)

Rv(i)=omega(i)*R*delta/upsilon;

Ra=2*nu*delta/upsilon;

lambda(i)=n0*Ra^m0*(1+(Rv(i)/Ra)^2)^((1+m0)/2);

sigma(i)=lambda(i)*l/delta;

B(i)=2-((Rv(i)/Ra)^2-m0)/((Rv(i)/Ra)^2+1);

E(i)=(1+z)/2/(1+z+2*sigma(i));

u0(i)=2*sigma(i)^2*E(i)*(1-m0)/(1+z+2*sigma(i));

u1(i)=2*sigma(i)^2*(E(i)/sigma(i)+B(i)*(1/6+E(i))/2)/(1+z+2*sigma(i));

u2(i)=sigma(i)*(1/6+E(i))/(1+z+2*sigma(i));

u3(i)=pi*R*p/lambda(i);

% 密封力Muszynska模型系數(shù)

% K=K0*(1-e^2)^(-n)

% D=D0*(1-e^2)^(-n)

% tau=tau0*(1-e)^b

% e=sqrt(x^2+y^2)/c

K0(i)=u0(i)*u3(i);

D0(i)=u1(i)*u3(i)*T;

mf(i)=u2(i)*u3(i)*(T^2);

M(i)=mf(i)+m;

% 系統(tǒng)結(jié)構(gòu)動力學(xué)方程中的參數(shù)

data=[m,g,c,n,b,mf(i),M(i),D0(i),De,tau0,Ke,K0(i),omega(i),e];

%初始條件

f0=[X0;0;Y0;0];

t0=0;h=0.001;tf=600;

tn=[t0:h:tf];

[t,f]=ode45(@rhs,tn,f0,[],data);

%求出的是無量綱化的位移、速度

%現(xiàn)轉(zhuǎn)化為實際的位移、速度

x=f(:,1);y=f(:,3);

dx=f(:,2);dy=f(:,4);

nn=(tf-t0)/h+1;nnn=300000;%nnn為除去瞬態(tài)響應(yīng)后的第一點數(shù)，nn為采集的總點數(shù)

nnnn=400000;

for i=nnn:nn

fprintf(fidd,'%6.4f??%12.8f\n',t(i),x(i));hold on;

end

fclose(fidd);

for j=nnn:nn

fprintf(fiddd,'%12.8f %12.8f\n',x(j),y(j));hold on;

end

fclose(fiddd);

% % 求龐加萊映射點

[A,jj]=max(x(nnn:nnnn));

jj=jj+nnn-1;

TT=2*pi;

for k=1:nn

tt(k)=(jj-1)*h+k*TT;

if tt(k)>tf

break;

end

for kk=nnn:nn

if tt(k)>(kk-1)*h&tt(k)

l1=(tt(k)-(kk+1)*h)*(tt(k)-kk*h)/(2*h^2);

l2=(tt(k)-(kk-1)*h)*(tt(k)-(kk+1)*h)/(-h^2);

l3=(tt(k)-(kk-1)*h)*(tt(k)-kk*h)/(2*h^2);

xx(k)=l1*x(kk-1)+l2*x(kk)+l3*x(kk+1);

dx(k)=l1*dx(kk-1)+l2*dx(kk)+l3*dx(kk+1);

fprintf(fid,'%12.8f %12.8f\n',xx(k),dx(k));hold on;

end

end

end

end

fclose(fid);

總結(jié)

以上是生活随笔為你收集整理的matlab 混沌分叉图,如何做出poincare 图或混沌分叉图的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。