计算物理 混沌 单摆.docx
- 文档编号:5483201
- 上传时间:2022-12-17
- 格式:DOCX
- 页数:10
- 大小:17KB
计算物理 混沌 单摆.docx
《计算物理 混沌 单摆.docx》由会员分享,可在线阅读,更多相关《计算物理 混沌 单摆.docx(10页珍藏版)》请在冰豆网上搜索。
计算物理混沌单摆
%例1:
单摆的不同周期的位移曲线
functiondjdb
[t1,w1]=ode45(@f,[0:
3*pi/200:
6*pi],[pi/7,0],[]);
[t2,w2]=ode45(@f,[0:
3*pi/200:
6*pi],[pi/3,0],[]);
plot(t1,w1(:
1),'-',t2,w2(:
1),'k:
');
xlabel('t');
ylabel('\theta');
%-----------------------------
functionydot=f(t,y)
ydot=[y
(2);
-sin(y
(1))];
%-----------------------------
%例2:
单摆的周期与摆角的关系
functiondjdb
theta=linspace(0.02,pi-0.02,40);
T=[];
options=odeset('Events',@events);%开启事件判断功能
%解不同的初始角度下的周期值
fork=1:
40;
[t,u]=ode45(@f,[0:
40],[theta(k),0],options);
T=[T,2*t(end)];
end
%解析解
k2=sin(theta./2).^2;
[K,E]=ellipke(k2);
plot(theta,T,theta,4*K,'*')
title('周期与摆角的关系');
xlabel('摆角');ylabel('周期');
%-----------------------------
functionydot=f(t,y)
ydot=[y
(2);
-sin(y
(1))];
%-----------------------------
%判断事件的变量为y
(2)=0,设置为y
(2)从负值增加到零时,中断计算
function[value,isterminal,direction]=events(t,y)
value=y
(2);
isterminal=1;direction=1;
%例3:
无阻尼无驱动时单摆的相图
%简单的程序
ezplot('0.5*y^2-cos(x)+0.8')
holdon
ezplot('0.5*y^2-cos(x)+0')
ezplot('0.5*y^2-cos(x)-0.6')
ezplot('0.5*y^2-cos(x)-1')
ezplot('0.5*y^2-cos(x)-1.4')
%复杂些的程序
%标注文字
plot([4.5,5.2],[0.8,0.8],'g',[4.5,5.2],[0,0],'r',[4.5,5.2],[-0.8,-0.8],'b');
text(5.3,0.8,'E<2mgl');
text(5.3,0,'E=2mgl');
text(5.3,-0.8,'E>2mgl');
xlabel('θ');ylabel('dθ/dt');
holdon
%能量方程
ydot=inline('sqrt(abs(E-1+cos(x)))','x','E');
e=[3,2.5,2,1.5,1,0.5,0.3,0.1];
%不同能量下的相图
fork=1:
8
ifk>3%对应E<2mgl
Q{k}=acos(1-e(k));
X=linspace(-Q{k},Q{k},300);
y=ydot(X,e(k));
plot(X,y,'g',X,-y,'g')
elseifk==3%对应E=2mgl
X=linspace(-2*pi,2*pi,300);
y=ydot(X,e(k));
plot(X,y,'r',X,-y,'r')
else%对应e>2mgl
X=linspace(-2*pi,2*pi,300);
y=ydot(X,e(k));
plot(X,y,'b',X,-y,'b')
end
end
holdoff
%例4:
有阻尼有驱动的三维相图
functiondbyd
globalafu
u=2/3;a=0.25;ZQ=3*pi;f=0.8;
[T,Y]=ode45(@dby,[0:
ZQ/100:
10*ZQ],[-0.8,2,u]);
%画相图
plot3(Y(:
1),Y(:
2),Y(:
3))
view(-95,60)
functionydot=dby(t,y)
globalafu
ydot=[y
(2);
-sin(y
(1))-2*a*y
(2)+f*cos(y(3))
u];
%例5:
对称性破缺、倍周期分岔与混沌
functiondbyd
globalafu
u=2/3;a=0.5;ZQ=3*pi;f=1.098;
[T,Y]=ode45(@dby,[0:
ZQ/100:
30*ZQ],[-0.8,2]);
subplot(2,1,1)%画相图
plot(Y(2500:
end,1),Y(2500:
end,2))
subplot(2,1,2)%画庞加莱截面图
xx=[];yy=[];
%forj=2500:
100:
length(Y(1:
end,1))-1
forj=2500:
100:
3000
xx=[xx,Y(j,1)];yy=[yy,Y(j,2)];
end
plot(xx,yy,'r.')
functionydot=dby(t,y)
globalafu
ydot=[y
(2);
-sin(y
(1))-a*y
(2)+f*cos(u*t)];
%例6:
分岔图——速度与f的关系
functiondbyd
globalafu
u=2/3;a=0.5;ZQ=3*pi;ff=1.05:
0.00005:
1.095;
FF=[];
YY=[];
fork=1:
length(ff)
f=ff(k);
[T,Y]=ode23(@dby,[0:
ZQ/100:
30*ZQ],[-0.8,2]);
FF=[FF,f];
YY=[YY;Y(2501:
100:
end,2)'];
end
plot(FF,YY,'r.','markersize',1)
functionydot=dby(t,y)
globalafu
ydot=[y
(2);
-sin(y
(1))-a*y
(2)+f*cos(u*t)];
%例7:
无阻尼无驱动时倒摆运动的相图
ezplot('y^2+x^4/2-x^2',[-1.6,1.6])
holdon
ezplot('y^2+x^4/2-x^2+0.2')
ezplot('y^2+x^4/2-x^2+0.4')
ezplot('y^2+x^4/2-x^2-0.4')
%例8:
倒摆运动的分岔图
functiondbd
globald
x0=0.1;v0=0.1;
d0=0.5:
0.002:
1.5;
axis([0.51.501.5])
holdon
forj=1:
length(d0)
d=d0(j);
[t,u]=ode45(@dbdfun,[0:
2*pi/60:
60*pi],[x0,v0]);
plot(d,u(901:
60:
1800,2),'r.');
end
functionydot=dbdfun(t,y)
globald
r=1;w=1;
ydot=[y
(2);
-y
(1)^3+y
(1)-d*y
(2)+r*cos(w*t)];
%例9:
倒摆的混沌运动及各种研究方法
functiondb
globald
%设v0有微小的变化,比较解的变化情况
x0=0.1;v0=0.1;d=0.78;
[t,u]=ode45(@dbfun,[0:
0.01:
100],[x0,v0]);
[t1,u1]=ode45(@dbfun,[0:
0.01:
100],[x0,v0-0.001]);
figure
plot(t,u(:
1),'r',t1,u1(:
1),'g')
xlabel('时间');
ylabel('摆角');
title('混沌状态下初条件有微小差异会形成的两条分开的曲线');
%当d=1.5,为周期1吸引子;当d=1.35为周期2吸引子;当d=1.15为奇怪吸引子
%可以改变d值,以观察不同的情况
d0=[1.5,1.35,1.15,1.05];
str{1}='庞加莱截面―周期1吸引子';
str{2}='庞加莱截面―周期2吸引子';
str{3}='庞加莱截面―奇怪吸引子';
str{4}='庞加莱截面―周期3吸引子';
forj=1:
4
d=d0(j);
[t,u]=ode45(@dbfun,[0:
2*pi/300:
200*pi],[x0,v0]);
figure
set(gcf,'unit','normalized','Position',[0.040.040.940.8]);
subplot(2,2,1)%位移曲线
plot(t,u(:
1))
title('位移曲线');
axis([0,150,-2.5,2.5]);
xlabel('t');ylabel('x');
subplot(2,2,2)%相图(奇怪吸引子)
plot(u(20000:
end,1),u(20000:
end,2))
title('相图');
axis([-22-1.51.5])
xlabel('x');ylabel('v');
Y=fft(u(:
1));%傅里叶功率分析
Y
(1)=[];n=length(Y);m=fix(n/2);%fix向零取整
power=abs(Y(1:
m)).^2/n^2;%功率,振幅的平方
freq=300*(1:
n/2)./(n*2*pi);%频率
subplot(2,3,4)
plot(freq,power)
axis([00.600.15])
title('功率谱');
xlabel('频率/Hz');ylabel('功率/w');
subplot(2,3,5)%庞加莱截面
plot(u(2000:
300:
30000,1),u(2000:
300:
30000,2),'r.');
axis([-22-1.51.5])
title(str{j});
subplot(2,3,6)%实物模拟图
h=plot([0,sin(x0)],[0,cos(x0)],'o-','erasemode','xor');
axis([-11-11])
title('倒摆运动模拟');
fori=25000:
30000
set(h,'xData',[0,sin(u(i,1))],'yData',[0,cos(u(i,1))]);
drawnow
end
end
%---------------------
functionydot=dbfun(t,y)
globald
r=1;w=1;
ydot=[y
(2);
-y
(1)^3+y
(1)-d*y
(2)+r*cos(w*t)];
%例10:
强迫VDP方程的各种研究方法
functionzjzd
globaluuwv
w=0.44;T=2*pi/w;v=1;
u=[0.851.030.61.0732];
str{1}='庞加莱截面―周期1吸引子';
str{2}='庞加莱截面―周期2吸引子';
str{3}='庞加莱截面―不变环面吸引子';
str{4}='庞加莱截面―奇怪吸引子';
forj=1:
4
figure
uu=u(j);
[t,y]=ode45(@vdpfun,[0:
T/1000:
100*T],[4,4]);
set(gcf,'unit','normalized','Position',[0.040.040.940.8]);
subplot(2,2,1)%位移曲线
plot(t(90000:
end),y(90000:
end,1))
title('位移曲线');
xlabel('x');ylabel('v');
subplot(2,2,2)%相图(奇怪吸引子)
plot(y(3000:
end,1),y(3000:
end,2))
title('相图');
axis([-33-44])
xlabel('x');ylabel('v');
Y=fft(y(:
1));%傅里叶功率分析
Y
(1)=[];n=length(Y);m=fix(n/2);%fix向零取整
power=abs(Y(1:
m)).^2/n^2;%功率,振幅的平方
dt=t
(2)-t
(1);
freq=1/dt*(1:
n/2)./n;%正确计算频率
subplot(2,2,3)
plot(freq,power)
axis([00.400.06])
title('功率谱');
xlabel('频率/Hz');ylabel('功率/w');
subplot(2,2,4)%庞加莱截面
fori=5000:
1000:
100000
plot(y(i,1),y(i,2),'r.');
holdon
end
axis([-31-11])
title(str{j});
end
%---------------------
functionydot=vdpfun(t,y)
globaluuwv
x0=1;w0=1;
ydot=[y
(2);
uu*(x0^2-y
(1)^2)*y
(2)-y
(1)*w0^2-v*cos(w*t)];
%例11:
洛伦兹方程的相图
functionlor
[t,y]=ode45(@lorfun,[0,31],[0.20,0.1,eps]);
subplot(2,2,1)%
plot3(y(:
1),y(:
2),y(:
3),y(1,1),y(1,2),y(1,3),'*')
xlabel('x');ylabel('y');zlabel('z');
view(3)
subplot(2,2,2)%
plot(y(:
1),y(:
2),y(1,1),y(1,2),'*')
xlabel('x');ylabel('y');
subplot(2,2,3)
plot(y(:
1),y(:
3),y(1,1),y(1,3),'*')
xlabel('x');ylabel('z');
subplot(2,2,4)%
plot(y(:
2),y(:
3),y(1,2),y(1,3),'*')
xlabel('y');ylabel('z');
%---------------------
functionydot=lorfun(t,y)
r=28;
ydot=[-10*y
(1)+10*y
(2);
r*y
(1)-y
(1)*y(3)-y
(2);
y
(1)*y
(2)-8/3*y(3)];
%例12:
洛伦兹方程的截面图
functionlor
opts=odeset('Events',@events);
[t,y,te,ye,ie]=ode45(@lorfun,[0,600],[1,0,0],opts);
plot(ye(:
2),ye(:
3),'.')
xlabel('y');ylabel('z');
%---------------------
functionydot=lorfun(t,y)
r=25;
ydot=[-10*y
(1)+10*y
(2);
r*y
(1)-y
(1)*y(3)-y
(2);
y
(1)*y
(2)-8/3*y(3)];
function[value,isterminal,direction]=events(t,y)
value=y
(1);isterminal=0;direction=0;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算物理 混沌 单摆 计算 物理
![提示](https://static.bdocx.com/images/bang_tan.gif)