计算物理 混沌 单摆.docVIP

%例1：单摆的不同周期的位移曲线 function djdb [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); %----------------------------- function ydot=f(t,y) ydot=[y(2); -sin(y(1))]; %----------------------------- %例2：单摆的周期与摆角的关系 function djdb theta=linspace(0.02,pi-0.02,40); T=[]; options=odeset(Events,@events);%开启事件判断功能 %解不同的初始角度下的周期值 for k=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(周期); %----------------------------- function ydot=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) hold on 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,E2mgl); text(5.3,0,E=2mgl); text(5.3,-0.8,E2mgl); xlabel(θ); ylabel(dθ/dt); hold on %能量方程 ydot=inline(sqrt(abs(E-1+cos(x))),x,E); e=[3, 2.5, 2, 1.5,1, 0.5, 0.3, 0.1]; %不同能量下的相图 for k=1:8 if k3 %对应E2mgl 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) elseif k==3 %对应E=2mgl X=linspace(-2*pi,2*pi,300); y=ydot(X,e(k)); plot(X,y,r,X,-y,r) else %对应e2mgl X=linspace(-2*pi,2*pi,300); y=ydot(X,e(k)); plot(X,y,b,X,-y,b) end end hold off %例4：有阻尼有驱动的三维相图 function dbyd global a f u 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) function ydot=dby(t,y) global a f u ydot=[y(2); -sin(y(1))-2*a*y(2) + f*cos(y(3)) u]; %例5：对称性破缺、倍周期分岔与混沌 function dbyd global a f u 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) %画相图 pl