单摆运动的描述（1）无阻尼单摆（小角度）20+*sin()0θωθ= 上式中令 sin()θθ=，201ω=得到如下方程：+0θθ=上述方程即为相图的方程，可由此方程画出无阻尼单摆在小角度下的相图：代码如下：%w0=2%E=2时syms x y ;%x 表示角度,y 表示角速度ezplot('x.^2+4*y.^2-4'),hold on%E=3时syms x y ;ezplot('x.^2+4*y.^2-6'),hold on%E=4时syms x y ;ezplot('x.^2+4*y.^2-8'),hold on%E=0.5时syms x yezplot('x.^2+4*y.^2-1'),hold onxlabel('角度')ylabel('角速度')title('无阻尼小角度单摆运动相图')上图中不同的同心椭圆表示在不同的能量下单摆的运动相图，在画上图时，令02ω=，改变能量E 得到一簇同心椭圆。

改变0ω会改变椭圆的形状，当01ω=时，椭圆变成圆。

下面时无阻尼小角度单摆的运动轨迹分析：此时只要求解上述的微分方程，然后改变其中的初始条件00(,)θω即可，其中求解微分方程的代码如下：%w0=1时%初始角度为pi/4时dsolve('D2y+y=0','y(0)=pi/4,Dy(0)=0','t')%用y 表示角度，Dy 表示角速度 %初始角度为pi/3时dsolve('D2y1+y1=0','y1(0)=pi/3,Dy1(0)=0','t')%此时令y1为角度%初始角度为pi/2时dsolve('D2y2+y2=0','y2(0)=pi/2,Dy2(0)=0','t')%此时用y3表示角度 画图的代码如下：%初始角度为pi/4时t=0:pi/50:4*pi;y=(pi*cos(t))/4;plot(t,y),hold on%初始角度为pi/3时y=(pi*cos(t))/3;plot(t,y,'r'),hold on%初始角度为pi/2时y=(pi*cos(t))/2;plot(t,y,'g'),hold onxlabel('时间')ylabel('角度')title('无阻尼小角度单摆在不同初始角度下的运动轨迹')legend('初始角度为pi/4的图','初始角度为pi/3的图','初始角度为pi/2的图') 出的图如下：当初始角度固定不变，改变初始角速度时，也会画出图，此时将初始角度固定为pi/4，画这些图首先需要求解运动方程，求解方法和上述相同，求解微分方程的代码如下：%研究初始角度相同，初始角速度不同的时候单摆的运动轨迹%此时都令初始的角度为pi/4%初始角速度为0时dsolve('D2y+y=0','y(0)=pi/4,Dy(0)=0','t')%用y表示角度，Dy表示角速度%初始角速度为1时dsolve('D2y+y=0','y(0)=pi/4,Dy(0)=1','t')%用y表示角度，Dy表示角速度%初始角速度为2时dsolve('D2y+y=0','y(0)=pi/4,Dy(0)=2','t')%用y表示角度，Dy表示角速度%初始角速度为3时dsolve('D2y+y=0','y(0)=pi/4,Dy(0)=3','t')%用y表示角度，Dy表示角速度下面是画图的代码：t=0:pi/50:4*pi;%w0=0y=(pi*cos(t))/4;plot(t,y),hold on%w0=1y=sin(t) + (pi*cos(t))/4;plot(t,y,'r'),hold on%w0=2y=2*sin(t) + (pi*cos(t))/4;plot(t,y,'g'),hold on%w0=3y=3*sin(t) + (pi*cos(t))/4;plot(t,y,'c'),hold onxlabel('时间')ylabel('角度')title('相同初始角度，不同初始角速度下的运动轨迹')legend('初始角速度为0的图','初始角速度为1的图','初始角速度为2的图','初始角速度为3的图')画出的图如下：（2）倒立摆分析 倒立摆的方程为-=0ϕϕ ，和小角度摆相同的是，同样通过求解方程得出运动轨迹，根据上式直接画出相图。

画图的代码如下：%倒立摆的运动相图绘制%w0=1syms y x ;%E=0时h=ezplot('y^2-x^2'),hold onset(h,'color','c')%E=1时ezplot('y^2-x^2-2'),hold on%E=2时h=ezplot('y^2-x^2-4'),hold onset(h,'color','r')%E=3时h=ezplot('y^2-x^2-6'),hold onset(h,'color','b')%E=4时h=ezplot('y^2-x^2-8'),hold onset(h,'color','k')xlabel('角度')ylabel('角速度')title('倒立摆的相图')legend('E=0的图','E=1的图','E=2的图','E=3的图','E=4的图')画出的图如下：其中E=0的图表示渐近线。

下面分析倒立摆的运动轨迹：<1>固定初始角速度为0，改变初始角度求解微分方程代码如下：%倒立摆的运动轨迹求解%初始角速度为0,改变初始角度%初始角度为pi/4dsolve('D2y-y=0','y(0)=pi/4,Dy(0)=0','t')%初始角度为pi/3dsolve('D2y-y=0','y(0)=pi/3,Dy(0)=0','t')%初始角度为pi/6dsolve('D2y-y=0','y(0)=pi/6,Dy(0)=0','t')画图的代码如下：%初始角速度为0%初始角度为pi/4t=0:0.01:2;y=(pi*exp(t))/8 + (pi*exp(-t))/8;plot(t,y),hold on%初始角度为pi/3y=(pi*exp(t))/6 + (pi*exp(-t))/6;plot(t,y,'r'),hold on%初始角度为pi/6y=(pi*exp(t))/12 + (pi*exp(-t))/12;plot(t,y,'g')xlabel('时间')ylabel('角度')title('倒立摆在不同初始角度下的运动轨迹')legend('初始角度为pi/4的图','初始角度为pi/3的图','初始角度为pi/6的图')画出的图如下：<2>固定初始角度为pi/4,改变初始角速度求解微分方程代码如下：%初始角度为pi/4,改变初始角速度%角速度为1dsolve('D2y-y=0','y(0)=pi/4,Dy(0)=1','t')%角速度为2dsolve('D2y-y=0','y(0)=pi/4,Dy(0)=2','t')%角速度为3dsolve('D2y-y=0','y(0)=pi/4,Dy(0)=3','t')画图的代码如下：%初始角度为pi/4,改变初始角速度%w0=1t=0:0.01:2;y=exp(t)*(pi/8 + 1/2) + exp(-t)*(pi/8 - 1/2);plot(t,y),hold on%w0=2y=exp(t)*(pi/8 + 1) + exp(-t)*(pi/8 - 1);plot(t,y,'r'),hold on%w0=3y=exp(t)*(pi/8 + 3/2) + exp(-t)*(pi/8 - 3/2);plot(t,y,'g')xlabel('时间')ylabel('角度')title('倒立摆在不同初始角速度下的运动轨迹')legend('角速度为1','角速度为2','角速度为3')画出的图如下：（3）有阻尼单摆有阻尼单摆的运动方程为：+2**+=0θβθθ其中 是阻尼系数。

有阻尼单摆分为三种情况，若阻尼，临界阻尼和过阻尼，下面画图来显示这三种情况的轨迹：画运动轨迹就是求解上述微分方程的问题，代码如下：%阻尼单摆的运动轨迹绘制%若阻尼下，阻尼系数为0.1dsolve('D2y+0.2*Dy+y','y(0)=pi/4,Dy(0)=0')t=0:0.01:60;y=(pi*exp(-t/10).*cos((3*11^(1/2)*t)/10))/4 +(11^(1/2)*pi*exp(-t/10).*sin((3*11^(1/2)*t)/10))/132plot(t,y),hold on%阻尼系数逐渐增大，到临界阻尼状态,阻尼系数为1dsolve('D2y+2*Dy+y','y(0)=pi/4,Dy(0)=0')y=(pi*exp(-t))/4 + (pi*t.*exp(-t))/4;plot(t,y,'r')%阻尼系数逐渐增大，到过阻尼状态,阻尼系数为2dsolve('D2y+4*Dy+y','y(0)=pi/4,Dy(0)=0')y=exp(t*(3^(1/2) - 2))*(pi/8 + (pi*3^(1/2))/12) +(3^(1/2)*pi*exp(-t*(3^(1/2) + 2))*(3^(1/2) - 2))/24;plot(t,y,'g')xlabel('时间')ylabel('角度')title('有阻尼单摆在不同阻尼系数下的运动轨迹')legend('弱阻尼下的运动轨迹','临界阻尼下的运动轨迹','过阻尼下的运动轨迹')画出的图如下：下面是相图绘制，相图就是绘制θ的一阶导和θ的关系式。