matlab中 n 039 是什么意思,Matlab讨论区

编程入门 行业动态 更新时间:2024-10-25 02:28:51

matlab中 n 039 是什么意思,Matlab<a href=https://www.elefans.com/category/jswz/34/1758357.html style=讨论区"/>

matlab中 n 039 是什么意思,Matlab讨论区

ww.jpg (86.06 KB, 下载次数: 0)

2013-8-25 12:00 上传

这是系统随参数u的分岔图

function dbf4

w=1.0;u=0.04;Pm=0.1;Pa=0.18;g=0.1;

T=2*pi/w;h=T/100;m=100;t0=0.0;ep=0.000001;

y0=[0;0];

options=odeset;options.RelTol=1e-6;options.AbsTol=1e-6;

ts=[0:T/100:250*T];

[t,dx]=ode45(@sb,ts,y0,options);%积分250周期作为迭代初值

x(1)=dx(end,1);x(2)=dx(end,2);a=1;b=1;

c(1)=h/2;c(2)=c(1);c(5)=c(1);

c(3)=h;c(4)=h;s(1)=x(1);s(2)=x(2);

v=1;

while v==1

t=t0;

x(1)=s(1);x(2)=s(2);

x1(1)=1.0;x1(2)=0.0;

x2(1)=0.0;x2(2)=1.0;

for i=1:4*m  %m表示迭代一周期的次数

t1=t;ts=t;

for ii=1:2

p(ii)=x(ii);w(ii)=x(ii);

end

for jj=1:4

f(1)=p(2);

if p(1)>=1

f(2)=-2*u*p(2)-(1+g*cos(t))*(p(1)-1)+Pm+Pa*cos(t);

elseif p(1)<1&p(1)>-1

f(2)=-2*u*p(2)++Pm+Pa*cos(t);

else p(1)<=-1

f(2)=-2*u*p(2)-(1+g*cos(t))*(p(1)+1)+Pm+Pa*cos(t);

end

for ii=1:2

p(ii)=c(jj).*f(ii)+w(ii);

x(ii)=c(jj+1).*f(ii)/3.0+x(ii);%在时间间隔h上利用RK算法积分一次

end

end

t=t1;ts=t;

for ii=1:2

p(ii)=x1(ii);

w(ii)=x1(ii);

end

for jj=1:4

f(1)=p(2);

if x(1)>=1

f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);

elseif x(1)<1&x(1)>-1

f(2)=-2*u*p(2);

else x(1)<=-1

f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);

end

t=ts+c(jj);

for ii=1:2

p(ii)=c(jj).*f(ii)+w(ii);

x1(ii)=c(jj+1).*f(ii)/3.0+x1(ii);%在时间间隔h上利用RK算法积分一次

end

end

t=t1;ts=t;

for ii=1:2

p(ii)=x2(ii);

w(ii)=x2(ii);

end

for jj=1:4

f(1)=p(2);

if x(1)>=1

f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);

elseif x(1)<1&x(1)>-1

f(2)=-2*u*p(2);

else x(1)<=-1

f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);

end

t=ts+c(jj);

for ii=1:2

p(ii)=c(jj).*f(ii)+w(ii);

x2(ii)=c(jj+1).*f(ii)/3.0+x2(ii);%在时间间隔h上利用RK算法积分一次

end

end

b=1;N(a,b)=t;

b=b+1;N(a,b)=x(1);

b=b+1;N(a,b)=x(2);

a=a+1;

N

end

r(1)=s(1)-x(1);

r(2)=s(2)-x(2);

dr(1,1)=1.0-x1(1);

dr(1,2)=-x2(1);

dr(2,1)=-x1(2);

dr(2,2)=1.0-x2(2);

I=[1 0;0 1];

P=I-dr;%雅可比矩阵

D=eig(P);%特征乘子

d=dr(1,1).*dr(2,2)-dr(1,2).*dr(2,1);

b=-r(2).*dr(1,2)+r(1).*dr(2,2);

e=-r(1).*dr(2,1)+r(2).*dr(1,1);

ds(1)=b/d;

ds(2)=e/d;

if abs(r(1))

break

else

s(1)=s(1)-ds(1);

s(2)=s(2)-ds(2);

end

end

D

function dx=sb(t,x)

w=1.0;u=0.04;Pm=0.1;Pa=0.18;g=0.1;

if x(1)>=1;

dx=[x(2);-2*u*x(2)-(1+g*cos(w*t))*(x(1)-1)+Pm+Pa*cos(w*t)];

elseif x(1)<1&x(1)>-1;

dx=[x(2);-2*u*x(2)+Pm+Pa*cos(w*t)];

else x(1)<=-1;

dx=[x(2);-2*u*x(2)-(1+g*cos(w*t))*(x(1)+1)+Pm+Pa*cos(w*t)];

end

这是齿轮系统打靶法程序(4周期打靶)

按u从大到小变化

u=0.048时系统开始变为4周期运动

u=0.034时系统开始变为8周期运动

程序结果为:这期间特征乘子都在单位圆内,表示为恒定为4周期运动

与分岔图不符(分岔图明显有分岔行为);还有就是这期间随u变化,特征乘子还出现共轭复数不止一次,比如:u=0.039时,特征乘子为-0.2197+0.3042i与-0.2197-0.3042i;u=0.038时,特征乘子为:-0.2129+0.2864i与-0.2129-0.2864i;

不知道出现共轭复数正不正常?或是程序有问题,请各位指正!

更多推荐

matlab中 n 039 是什么意思,Matlab讨论区

本文发布于:2024-02-11 19:18:34,感谢您对本站的认可!
本文链接:https://www.elefans.com/category/jswz/34/1682877.html
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。
本文标签:讨论区   matlab   Matlab

发布评论

评论列表 (有 0 条评论)
草根站长

>www.elefans.com

编程频道|电子爱好者 - 技术资讯及电子产品介绍!