用matlab求解单摆模型

【实例简介】
在计算物理运用matlab软件进行编程，从而有效处理单摆模型，得到从周期运动到混沌图像
椭圆点:坐标原点(6=0,d0t=0),对 分界线 应K=V=E0,单摆的静止平衡点,附近 轨线为圆形或椭圆形的闭合轨道,轨道上 各点能量相等(等能轨道) 双曲奇点(鞍点):单摆倒置平衡点(O= vO 士丌,O6/t=0),附近相轨线为双曲线 分界线:从0=-丌,d/dt=0到0= 丌,O/Ot=0或相反的连线 十丌 E<max(V).轨线在分界线内,是闭合回线,单摆作周期振动 E≥max(V),轨线在分界线以外,轨道不闭合,单摆作旋转运动 3柱面上的单摆相轨线 相图横坐标是以2丌为周期的, 摆角土丌是单摆的同一个倒立位 双曲点 6-± 椭圆点 置,把相图上G点与G点重迭 起时,就把相平面卷缩成一个柱 面。所有相轨线都将呈现在柱面 上 2.有阻尼无驱动情况 此时≠0,f=0,方程是 explor d20 de +23,+sin6=0 有阻尼时,由于能量被消耗,单摆振动 的振幅会逐渐变小 有阻尼单摆的相图 能量耗散使相轨线矢径对数衰 减。无论从那点出发,经若干 次旋转后趋向坐标原点,原点称 为吸引子,它把相空间的点吸引 过来,原点又称不动点。 任意振幅下的相图 1.整个相平面被通过鞍点G与G的轨线分成三个区域 2.在坐标原点附近轨线由于单摆振幅不断减小而形成向内旋转 3鞍点的位置仍在G与G处 运动单摆从倒立开始运动后,由于 能量耗散达不到原有高度。 +1( 轨线从一个鞍点出发到不了另一鞍 点,分界线被破坏了 十元 相流所有中间区域的相点流向坐标 原点。原点是该区域的不动点,是该 区域吸引子。左右两个区域也有相应 的吸引子,它们分别处在该图左(-2丌 多多 )和右(+2丌)两侧 3.有阻尼有驱动情况 此时有≠0.,f≠0下面的计算取β=1/4,=2/3),方程是 20 2+2B*4+sin 0-f cos wt 为了构成单摆的三维相空间,将它化成三个方程: do dt dt dt sin b+f cos wt 用三个变量0,p,9组成三维相空间,相 角p有周期性,把2mm和2(m+1)丌平面连接 起来,相空间形成一个轮胎面。原来园形 2+1元y Cn 轨线成了在轮胎面上的环绕线 在轮胎面上作截面(庞加莱截面),轨线 穿过它会留下一个点,它对应取某个常 位相值时,在θ,平面上的相点 3.1极限环吸引子 当驱动力较小时,如f=0.8,从不同的 初始条件出发,经过长时间运动最终都 落到同一个椭圆上。称之为极限环 左图初始条件是0=-2,O/t=2, 1.4 轨线由外向内旋转到极限环 右图初始条件是θ de dl 4 1.74 0.2,轨线由内向外旋转到极限环 3.2对称性破缺 当∫=1.03,从不同的初条件出发,得 到不同的蛋形的吸引子,它们左右反 射对称,原来的左右对称性被破坏。4 左图初始条件是0=-0.1,0/t=2, 1登 在庞加莱截面仍然是一个点 右图初始条件是O=-0.8,d0t=2 在庞加菜截面仍然是一个点 161B1别U 33倍周期分岔与混沌 0 age 2 94 155 1 19 1 161 晶边2自 f=1.65 f=1.082, =1.088 倍周期分岔, 4周期分岔 出现混沌 计算中取0=-0.8,l0/=2。 34相图与庞加莱截面程序 这个程序主要学习相图与截面图的画法 function byd LegATO f u=2/3;a=0.5;因Q=3*pi;f=1.089; [T,Y=ode45(@dby,[0:zQ/200:500*2Q,[-0.8,2]); figure p1ot(Y(31000:end,1),Y(31000:end,2)) figure forj=40001:200:1 ength(Y(1:end,1)-1 xX=LXX, Y(j, 1)1 yy=[yy,Y(j,2)]; end p1ot(x,yy,’,r2) function dot=dby (t, y) global a f ydot=Ly(2) sin(y(1))-a*y(2)+f*Co8(u*)]; 35大摆角单摆程序 学习事件开关函数的用法 function diddb figure axiS([-88-22]) hold on %标注文字 p1ot([4.5,5.2],[0.8,0.8],’g2,[4.5,5.2],[0,O], r),[4.5,5.2],[0.8~-0.8],b2) 七ext(5.3,0.8,E<2mg1); text(5.3,0,E=2mg13); text(53,-0.8,2E>2mg1) xlabel(0); ylabel('d edt,) %能量方程 ydot=inline('sqrt(abs(E-1+Cos(x)))2,'x, 'E,) e=[3,2.5,2,1.5,1,0.5,0.3,0.1]; %不同能量下的相图 for k=1: 8 ifk>3%对应E<2mg Qlk]=acos(1-e(k)) X=linspace(-Qtk, Qtk], 300) y=dot(x, e(k)) plot(X,y,'g,X,y,'g) e] seif k==3%对应E=2mg1 X=linspace(-2*pi, 2*pi, 300) y=dot(x,e(k)) plot(X,y,r,,X,y, ',) else %对应e>2mg1 X=linspace(-2*pi, 2*pi, 300) y=dot(x, e(k)) p1ot(x,y,b,X,-y,b’) end end hold off %解不同初始角度下的微分方程 [七1,1]=ode45(f,[0:0.001:6],[pi/7,0],[); [t2,w2]=ode45(@f,[0:0.001:6],[pi/3,0],[]) %画不同角度下的波形图 figure plot(七1,W1(:,1),t2,W2(:,1)); x1abe1(时间);y1abe1(摆角’); legend()小摆角’,大摆角); %画周期与摆角的关系 theta=linspace(pi/360, pi-01, 40) options= odeset( Events, Events);%开启事件判断功能 %解不同的初始角度下的周期值 fori=1:40; [,u]=ode45(@f,[O:0.001:20],[hea(i),0], options); T=[T,2*t(end)]; end figure plot(theta, T) 七i1e()周期与摆角的关系); x]abe1(摆角’);y1abe1(周期); function dot=f(t, y) dot=Ly(2); 9.8*sin(y(1))]; function [value, isterminal, direction]=events(t, y) value=y(2) isterminal=1: direction=1 第五节倒摆与杜芬方程 1.倒摆实 11倒摆实验演示 1.2倒摆的简化模型与运动方程 倒摆可以简化成右图的模型,它的运动 以用杜芬方程描述 de dx +h +r cos wt 改变运动阻尼,可以演示运动状态从 周期解到混沌的变化 3.杜芬 Duffing)方程 下面用波形图,相图,频谱图和庞加莱截面图(map图)研究系统的运 动 3.1无阻尼无驱动情形 x+x3=0 积分得 1/da 所以势能是 4 这时有三个平衡点,x=0(不稳定叶 平衡点)x=±1(是稳定平衡点) 【实例截图】
【核心代码】