利用matlab进行混沌摆仿真（双摆、三摆、多摆演示）

本文首发于 matlab爱好者 微信公众号，欢迎关注。

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

本文为学习非线性动力学时，搜索资料时发现的一个小例子，感觉很有趣，试着自己尝试一下。本来想用计算多体动力学的方法求解，但时间不足，智商有限，没有学会。只能用丑陋的公式法，直接进行求解，实属遗憾。在此挖个小坑，日后有机会参悟出计算法后，必定在这里更新。

0 前言

非线性线性存在于物理与工程中的各个领域。在机械中，存在着大量的非线性现象。这里通过双摆和三摆的例子，来展示一个混沌现象。

什么是混沌呢，用俗话说，差别很小的输入，输入进混沌系统后，过一段时间总能将两者区分开。再换句话说，混沌系统对误差非常敏感，随着时间的推移，计算结果的可信度将会越来越低。反正现在科普文章还是很多的，有的文章把各种玄学都搬出来了。但系统就在那里，大家盲人摸象，不同的人会摸到不同的地方，产生自己的感悟。我这里也不啰嗦搞这一套了。

好了，不啰嗦了前言了。（话说这么啰嗦的前言会不会直接劝退读者 (っ°Д°;)っ ）

本文的混沌摆模型都是二维平面内的运动的模型，也就是平面铰链连接的机构，不涉及三维运动以及球铰连接。本文的混沌摆没有考虑摩擦力。

1 拉格朗日方法建模

这里假设摆的长度和质量都一样。当然，不一样的话，也可以求，就是公式复杂一些，方法一样。

1.1 双摆模型

上图是一个随手画的双摆模型，由两个小棍组成。上面小棍顶端点被限制，小面小棍顶端点与上面小棍相连。现在只需要用两个变量就可以描述这个双摆模型的状态，一个是蓝色摆与垂线的夹角θ1，一个是黄色摆与垂线的夹角θ2。这两个摆的角速度，可以用$\frac{d\theta 1}{dt}$、$\frac{d\theta 2}{dt}$来表示（或者简写为$\dot{{\theta }_1}$和$\dot{{\theta }_2}$）。

之后根据这两个信息，我们就可以愉快的列方程了：

首先，每个小棍对质心求转动惯量：
$$I=\frac{1}{12}m{l}^2$$
小棍1的质心位置为：
$$\begin{gathered} x_1=\frac{l}{2}\sin {\theta }_1 \\ {y}_1=-\frac{l}{2}\cos {\theta }_1 \end{gathered}$$
小棍2的质心位置比较复杂，为：
$$\begin{gathered} x_2=l(\sin {\theta }_1+\frac{1}{2}\sin {\theta }_2) \\ {y}_2=-l(\cos {\theta }_1+\frac{1}{2}\cos {\theta }_2) \end{gathered}$$
可以看做是小棍1的末端坐标，再加上小棍2的偏移量叠加得到的。
知道了质心坐标后，我们就可以求出质心的运动速度：
小棍1的质心速度为：
$$\begin{gathered} d{x}_1=\frac{l}{2}\cos {\theta }_1\cdot \dot{{\theta }_1} \\ d{y}_1=\frac{l}{2}\sin {\theta }_1\cdot \dot{{\theta }_1} \end{gathered}$$
然后小棍2的质心速度为：
$$\begin{gathered} d{x}_2=l(\cos {\theta }_1\cdot \dot{{\theta }_1}+1/2\cdot \cos {\theta }_2\cdot \dot{{\theta }_2}) \\ d{x}_2=-l(\sin {\theta }_1\cdot \dot{{\theta }_1}+1/2\cdot \sin {\theta }_2\cdot \dot{{\theta }_2}) \end{gathered}$$
然后求拉格朗日量（详情见任意理论力学课本）：
$$\begin{gathered} L=T_{平动}+T_{转动}-U_{势能} \\ =\frac{1}{2}m(d{x}_1^2+d{y}_1^2)+\frac{1}{2}m(d{x}_2^2+d{y}_2^2)+\frac{1}{2}I{\dot{{\theta }_1}}^2+\frac{1}{2}I{\dot{{\theta }_2}}^2-mg(y_1+y_2) \end{gathered}$$

然后我们需要求解出角度与角动量p相关的方程组，方便后续求解。首先是角度变化方程，这个需要方程求解一下：
$$\begin{gathered} p_1=\frac{\partial L}{\partial \dot{{\theta }_1}} \\ {p}_2=\frac{\partial L}{\partial \dot{{\theta }_2}} \end{gathered}$$
得到：
$$\begin{gathered} \dot{{\theta }_1}=\frac{6(2p_1-3p_2\cos ({\theta }_1-{\theta }_2))}{m{l}^2\left(9{\cos }^2({\theta }_1-{\theta }_2)-16\right)} \\ \dot{{\theta }_2}=\frac{6(8p_2-3p_1\cos ({\theta }_1-{\theta }_2))}{m{l}^2\left(9{\cos }^2({\theta }_1-{\theta }_2)-16\right)} \end{gathered}$$
角动量变化方程比较简单，不需要联立求解：
$$\begin{gathered} \dot{p_1}=\frac{\partial L}{\partial {\theta }_1} \\ \dot{p_2}=\frac{\partial L}{\partial {\theta }_2} \end{gathered}$$
得到：
$$\begin{gathered} \dot{p_1}=\frac{1}{2}ml(3g\sin {\theta }_1+\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2)) \\ \dot{p_2}=\frac{1}{2}ml(-g\sin {\theta }_2+\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2)) \end{gathered}$$
然后就是用龙格库塔方法，求解下面的方程组就行了：
$$\{\begin{array}{c}\dot{{\theta }_1}=\frac{6(2p_1-3p_2\cos ({\theta }_1-{\theta }_2))}{m{l}^2\left(9{\cos }^2({\theta }_1-{\theta }_2)-16\right)}\\ \\ \dot{{\theta }_2}=\frac{6(8p_2-3p_1\cos ({\theta }_1-{\theta }_2))}{m{l}^2\left(9{\cos }^2({\theta }_1-{\theta }_2)-16\right)}\\ \\ \dot{p_1}=\frac{1}{2}ml(3g\sin {\theta }_1+\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2))\\ \\ \dot{p_2}=\frac{1}{2}ml(-g\sin {\theta }_2+\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2))\end{array}$$
当然以上过程，全部可以用matlab的符号求解功能，或者用mathematica进行求解，而且事实上强烈建议用软件进行求解，以免造成手算导致的失误。

1.2 三摆以及多摆模型

知道了双摆模型，三摆模型或者多摆模型就可以迎刃而解了。

如果采用程序法计算的话：
不同小棍质心点的坐标，用循环的方式可以有规律的生成。拉格朗日量也可以用循环的方式去构建。
之后求偏导以及解方程，也都可以自动化完成，非常的方便。

不过这里还是给出三摆的运动方程，方便各位读者校对：
各个小棍重心的位置方程如下，小棍1和小棍2同双摆：
$$\begin{gathered} x_1=\frac{l}{2}\sin {\theta }_1 \\ {y}_1=-\frac{l}{2}\cos {\theta }_1 \\ {x}_2=l(\sin {\theta }_1+\frac{1}{2}\sin {\theta }_2) \\ {y}_2=-l(\cos {\theta }_1+\frac{1}{2}\cos {\theta }_2) \\ {x}_3=l(\sin {\theta }_1+\sin {\theta }_2+\frac{1}{2}\sin {\theta }_3) \\ {y}_3=-l(\cos {\theta }_1+\cos {\theta }_2+\frac{1}{2}\cos {\theta }_3) \end{gathered}$$
小棍重心的变化为：
$$\begin{gathered} d{x}_1=\frac{l}{2}\cos {\theta }_1\cdot \dot{{\theta }_1} \\ d{y}_1=\frac{l}{2}\sin {\theta }_1\cdot \dot{{\theta }_1} \\ d{x}_2=l(\cos {\theta }_1\cdot \dot{{\theta }_1}+1/2\cdot \cos {\theta }_2\cdot \dot{{\theta }_2}) \\ d{x}_2=-l(\sin {\theta }_1\cdot \dot{{\theta }_1}+1/2\cdot \sin {\theta }_2\cdot \dot{{\theta }_2}) \\ d{x}_3=l(\cos {\theta }_1\cdot \dot{{\theta }_1}+\cos {\theta }_2\cdot \dot{{\theta }_2}+1/2\cdot \cos {\theta }_3\cdot \dot{{\theta }_3}) \\ d{x}_3=-l(\sin {\theta }_1\cdot \dot{{\theta }_1}+\sin {\theta }_2\cdot \dot{{\theta }_2}+1/2\cdot \sin {\theta }_3\cdot \dot{{\theta }_3}) \end{gathered}$$
这个三摆系统的拉格朗日量L为：
$$\begin{gathered} L=\frac{1}{2}m(d{x}_1^2+d{y}_1^2+d{x}_2^2+d{y}_2^2+d{x}_3^2+d{y}_3^2)+\frac{1}{2}I({\dot{{\theta }_1}}^2+{\dot{{\theta }_2}}^2+{\dot{{\theta }_3}}^2)-mg(y_1+y_2+y_3) \\ =\frac{1}{6}ml\left(9\dot{{\theta }_2}\dot{{\theta }_1}l\cos ({\theta }_1-{\theta }_2)+3\dot{{\theta }_3}\dot{{\theta }_1}l\cos ({\theta }_1-{\theta }_3)+3\dot{{\theta }_2}\dot{{\theta }_3}l\cos ({\theta }_2-{\theta }_3)+7{\dot{{\theta }_1}}^{2}l+4{\dot{{\theta }_2}}^{2}l+{\dot{{\theta }_3}}^{2}l+15g\cos ({\theta }_1)+9g\cos ({\theta }_2)+3g\cos ({\theta }_3)\right) \end{gathered}$$
同样，我们联立求解出角度的变化
$$\begin{gathered} p_1=\frac{\partial L}{\partial \dot{{\theta }_1}} \\ {p}_2=\frac{\partial L}{\partial \dot{{\theta }_2}} \\ {p}_3=\frac{\partial L}{\partial \dot{{\theta }_3}} \end{gathered}$$
得到
$$\begin{gathered} \dot{{\theta }_1}=\frac{6(9p_1\cos (2({\theta }_2-{\theta }_3))+27p_2\cos ({\theta }_1-{\theta }_2)-9p_2\cos ({\theta }_1+{\theta }_2-2{\theta }_3)+21p_3\cos ({\theta }_1-{\theta }_3)-27p_3\cos ({\theta }_1-2{\theta }_2+{\theta }_3)-23p_1)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)} \\ \dot{{\theta }_2}=\frac{6(27p_1\cos ({\theta }_1-{\theta }_2)-9p_1\cos ({\theta }_1+{\theta }_2-2{\theta }_3)+9p_2\cos (2({\theta }_1-{\theta }_3))-27p_3\cos (2{\theta }_1-{\theta }_2-{\theta }_3)+57p_3\cos ({\theta }_2-{\theta }_3)-47p_2)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)} \\ \dot{{\theta }_3}=\frac{6(21p_1\cos ({\theta }_1-{\theta }_3)-27p_1\cos ({\theta }_1-2{\theta }_2+{\theta }_3)-27p_2\cos (2{\theta }_1-{\theta }_2-{\theta }_3)+57p_2\cos ({\theta }_2-{\theta }_3)+81p_3\cos (2({\theta }_1-{\theta }_2))-143p_3)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)} \end{gathered}$$
另外三个角动量方程为：
$$\begin{gathered} \dot{p_1}=\frac{\partial L}{\partial {\theta }_1} \\ \dot{p_2}=\frac{\partial L}{\partial {\theta }_2} \\ \dot{p_3}=\frac{\partial L}{\partial {\theta }_3} \end{gathered}$$
计算得到：
$$\begin{gathered} \dot{p_1}=-\frac{1}{2}ml(3\dot{{\theta }_2}\dot{{\theta }_1}l\sin ({\theta }_1-{\theta }_2)+\dot{{\theta }_1}\dot{{\theta }_3}l\sin ({\theta }_1-{\theta }_3)+5g\sin ({\theta }_1)) \\ \dot{p_1}=-\frac{1}{2}ml(-3\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2)+\dot{{\theta }_2}\dot{{\theta }_3}l\sin ({\theta }_2-{\theta }_3)+3g\sin ({\theta }_2)) \\ \dot{p_1}=-\frac{1}{2}ml(\dot{{\theta }_1}\dot{{\theta }_3}l\sin ({\theta }_1-{\theta }_3)+\dot{{\theta }_2}\dot{{\theta }_3}l\sin ({\theta }_2-{\theta }_3)-g\sin ({\theta }_3)) \end{gathered}$$
因此，我们得到了下面这个6维一阶常微分方程组：
$$\{\begin{array}{c}\dot{{\theta }_1}=\frac{6(9p_1\cos (2({\theta }_2-{\theta }_3))+27p_2\cos ({\theta }_1-{\theta }_2)-9p_2\cos ({\theta }_1+{\theta }_2-2{\theta }_3)+21p_3\cos ({\theta }_1-{\theta }_3)-27p_3\cos ({\theta }_1-2{\theta }_2+{\theta }_3)-23p_1)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)}\\ \\ \dot{{\theta }_2}=\frac{6(27p_1\cos ({\theta }_1-{\theta }_2)-9p_1\cos ({\theta }_1+{\theta }_2-2{\theta }_3)+9p_2\cos (2({\theta }_1-{\theta }_3))-27p_3\cos (2{\theta }_1-{\theta }_2-{\theta }_3)+57p_3\cos ({\theta }_2-{\theta }_3)-47p_2)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)}\\ \\ \dot{{\theta }_3}=\frac{6(21p_1\cos ({\theta }_1-{\theta }_3)-27p_1\cos ({\theta }_1-2{\theta }_2+{\theta }_3)-27p_2\cos (2{\theta }_1-{\theta }_2-{\theta }_3)+57p_2\cos ({\theta }_2-{\theta }_3)+81p_3\cos (2({\theta }_1-{\theta }_2))-143p_3)}{m{l}^2(81\cos (2({\theta }_1-{\theta }_2))-9\cos (2({\theta }_1-{\theta }_3))+45\cos (2({\theta }_2-{\theta }_3))-169)}\\ \\ \dot{p_1}=-\frac{1}{2}ml(3\dot{{\theta }_2}\dot{{\theta }_1}l\sin ({\theta }_1-{\theta }_2)+\dot{{\theta }_1}\dot{{\theta }_3}l\sin ({\theta }_1-{\theta }_3)+5g\sin ({\theta }_1))\\ \\ \dot{p_1}=-\frac{1}{2}ml(-3\dot{{\theta }_1}\dot{{\theta }_2}l\sin ({\theta }_1-{\theta }_2)+\dot{{\theta }_2}\dot{{\theta }_3}l\sin ({\theta }_2-{\theta }_3)+3g\sin ({\theta }_2))\\ \\ \dot{p_1}=-\frac{1}{2}ml(\dot{{\theta }_1}\dot{{\theta }_3}l\sin ({\theta }_1-{\theta }_3)+\dot{{\theta }_2}\dot{{\theta }_3}l\sin ({\theta }_2-{\theta }_3)-g\sin ({\theta }_3))\end{array}$$

我这里提供一个简单的mathematica小程序，用来计算这个多摆方程用的（感觉matlab的兼容性应该会更好，毕竟之后我是用matlab求解，但是我个人不太习惯matlab的符号求解）：

下面是求解拉格朗日量的程序：

num = 3
II = 1/12*m*l^2;
Do[x[i] =
  Sum[l*Sin[th[k]], {k, 1, i}] - l/2*Sin[th[i]]
 , {i, 1, num}]
Do[y[i] =
  -Sum[l*Cos[th[k]], {k, 1, i}] + l/2*Cos[th[i]]
 , {i, 1, num}]
Do[dx[i] =
  Sum[l*Cos[th[k]]*dth[k], {k, 1, i}] - l/2*Cos[th[i]]*dth[i]
 , {i, 1, num}]
Do[dy[i] =
  Sum[l*Sin[th[k]]*dth[k], {k, 1, i}] - l/2*Sin[th[i]]*dth[i]
 , {i, 1, num}]
L = TrigReduce[
   1/2*m*(Sum[dx[i]^2, {i, 1, num}] + Sum[dy[i]^2, {i, 1, num}]) +
    1/2*II*(Sum[dth[i]^2, {i, 1, num}]) -
    m*g*(Sum[y[i], {i, 1, num}])];

然后求解三个角度

Solve[pth[1] == D[L, dth[1]] &&
  pth[2] == D[L, dth[2]] &&
  pth[3] == D[L, dth[3]], {dth[1], dth[2], dth[3]}]

最后计算三个动量

Do[dpth[i] = Simplify[D[L, th[i]]], {i, 1, num}]
dpth[1]
dpth[2]
dpth[3]

但是这个方程的复杂程度，感觉是和小棍的数量平方一个量级的（只是感觉，没有证明）。所以如果小棍的数量比较多，这种方法计算的方程量就会非常大。到那个时候，就不建议用理论计算出方程的方式去求解，建议用计算多体动力学之类的方式，进行数值求解。

2 多摆的数值求解

2.1 双摆

我这里没用ode45，用的是自己编的一个4阶龙格库塔方法进行求解的。

肉眼上看动图，看起来是没有什么问题的，运动连贯，混沌现象满足，能量、位置之类的也没有发散。

程序如下：

%双摆
clear
clc
close all
%输入
N=2;%双摆
m=1;
l=1;
g=9.8;
Input=[N,m,l,g];
%初始条件和时间设置
y0=[pi/2;pi/2;0;0];%这里全部是弧度值。分别代表[摆1与垂面夹角，摆2与垂面夹角，摆1角动量，摆2角动量]
h=1e-2;
x0=0:h:20;
%代入到ODE求解器中
[y1,Output]=ODE_RK4_hyh(x0,h,y0,Input);

%提取出角度
tN=size(y1,2);
th1=y1(1,:);
th2=y1(2,:);

%计算出关节坐标
CX1_A=zeros(1,tN);
CX1_B=CX1_A+l*sin(th1);
CY1_A=zeros(1,tN);
CY1_B=CY1_A-l*cos(th1);

CX2_A=CX1_B;
CX2_B=CX2_A+l*sin(th2);
CY2_A=CY1_B;
CY2_B=CY2_A-l*cos(th2);

%绘图
n=1;
figure()
set(gcf,'position',[488   342   400   300])
for k=1:4:length(x0) %这里4步一显示时间帧
    clf
    xlim([-2,2])
    ylim([-2,2])
    hold on
    %绘制摆
    plot([CX1_A(k),CX1_B(k)],[CY1_A(k),CY1_B(k)],'color','k','LineWidth',1.5)
    plot([CX2_A(k),CX2_B(k)],[CY2_A(k),CY2_B(k)],'color','k','LineWidth',1.5)
    %绘制轨线
    if k>200
        n=n+1;
    end
    Nm=k-n+1;
    %轨迹1
    F_color=[1,0,0];
    F_color=F_color*0.6+[1,1,1]*0.4*0.999;
    cdata=[linspace(1,F_color(1),Nm+1)',linspace(1,F_color(2),Nm+1)',linspace(1,F_color(3),Nm+1)'];
    cdata=reshape(cdata,Nm+1,1,3);
    if k>3
        patch([CX1_B(n:k),NaN],[CY1_B(n:k),NaN],1:Nm+1,'EdgeColor','interp','Marker','none',...
          'MarkerFaceColor','flat','CData',cdata,'LineWidth',1.5);
    end
    %轨迹2
    F_color=[0,0,1];
    F_color=F_color*0.6+[1,1,1]*0.4*0.999;
    cdata=[linspace(1,F_color(1),Nm+1)',linspace(1,F_color(2),Nm+1)',linspace(1,F_color(3),Nm+1)'];
    cdata=reshape(cdata,Nm+1,1,3);
    if k>3
        patch([CX2_B(n:k),NaN],[CY2_B(n:k),NaN],1:Nm+1,'EdgeColor','interp','Marker','none',...
          'MarkerFaceColor','flat','CData',cdata,'LineWidth',1.5);
    end
    hold off
    pause(0.05)
    %可以在这里添加输出动图的程序
end

function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn    ,yn       ,Input);
    K2=Fdydx(xn+h/2,yn+h/2*K1,Input);
    K3=Fdydx(xn+h/2,yn+h/2*K2,Input);
    K4=Fdydx(xn+h  ,yn+h*K3  ,Input);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
Output=[];
end

function dydx=Fdydx(x,y,Input)
%将原方程整理为dy/dx=F(y,x)的形式
%输入Input整理
m=Input(2);
l=Input(3);
g=Input(4);
%输入
th1=y(1);%角度1
th2=y(2);%角度2
pth1=y(3);%角动量1
pth2=y(4);%角动量2
%利用拉格朗日法得到的方程
M=l^2*m*(-16 + 9*cos(th1 - th2)^2);
dth1 = -6*(2*pth1 - 3*pth2*cos(th1 - th2))/M; 
dth2 = -6*(8*pth2 - 3*pth1*cos(th1 - th2))/M;
dpth1=-0.5*l*m*(3*g*sin(th1)+dth1*dth2*l*sin(th1-th2));
dpth2=0.5*l*m*(dth1*dth2*l*sin(th1-th2)-g*sin(th2));
%整理输出
dydx=zeros(4,1);
dydx(1)=dth1;
dydx(2)=dth2;
dydx(3)=dpth1;
dydx(4)=dpth2;
end

2.2 三摆

这是三摆模型，其实和双摆大同小异。

代码也和双摆差不多：

%三摆
clear
clc
close all
%输入
N=3;%双摆
m=1;
l=1;
g=9.8;
Input=[N,m,l,g];
%初始条件和时间设置
y0=[1/2*pi;1/2*pi;1/2*pi;0;0;0];%这里全部是弧度值。分别代表[摆1与垂面夹角，摆2与垂面夹角，摆1角动量，摆2角动量]
h=1e-2;
x0=0:h:20;
%代入到ODE求解器中
[y1,Output]=ODE_RK4_hyh(x0,h,y0,Input);

%提取出角度
tN=size(y1,2);
th1=y1(1,:);
th2=y1(2,:);
th3=y1(3,:);

%计算出关节坐标
CX1_A=zeros(1,tN);
CX1_B=CX1_A+l*sin(th1);
CY1_A=zeros(1,tN);
CY1_B=CY1_A-l*cos(th1);

CX2_A=CX1_B;
CX2_B=CX2_A+l*sin(th2);
CY2_A=CY1_B;
CY2_B=CY2_A-l*cos(th2);

CX3_A=CX2_B;
CX3_B=CX3_A+l*sin(th3);
CY3_A=CY2_B;
CY3_B=CY3_A-l*cos(th3);
%绘图
n=1;
figure()
set(gcf,'position',[488   342   400   300])
for k=1:4:1100
    clf
    xlim([-3,3])
    ylim([-3,3])
    hold on
    plot([CX1_A(k),CX1_B(k)],[CY1_A(k),CY1_B(k)],'color','k','linewidth',1)
    plot([CX2_A(k),CX2_B(k)],[CY2_A(k),CY2_B(k)],'color','k','linewidth',1)
    plot([CX3_A(k),CX3_B(k)],[CY3_A(k),CY3_B(k)],'color','k','linewidth',1)
    hold off
    
    %绘制轨线
    if k>200
        n=n+1;
    end
    Nm=k-n+1;
    %轨迹1
    F_color=[1,0,0];
    F_color=F_color*0.6+[1,1,1]*0.4*0.999;
    cdata=[linspace(1,F_color(1),Nm+1)',linspace(1,F_color(2),Nm+1)',linspace(1,F_color(3),Nm+1)'];
    cdata=reshape(cdata,Nm+1,1,3);
    if k>3
        patch([CX1_B(n:k),NaN],[CY1_B(n:k),NaN],1:Nm+1,'EdgeColor','interp','Marker','none',...
          'MarkerFaceColor','flat','CData',cdata,'LineWidth',1);
    end
    %轨迹2
    F_color=[0,0,1];
    F_color=F_color*0.6+[1,1,1]*0.4*0.999;
    cdata=[linspace(1,F_color(1),Nm+1)',linspace(1,F_color(2),Nm+1)',linspace(1,F_color(3),Nm+1)'];
    cdata=reshape(cdata,Nm+1,1,3);
    if k>3
        patch([CX2_B(n:k),NaN],[CY2_B(n:k),NaN],1:Nm+1,'EdgeColor','interp','Marker','none',...
          'MarkerFaceColor','flat','CData',cdata,'LineWidth',1);
    end
    %轨迹3
    F_color=[0,1,0];
    F_color=F_color*0.6+[1,1,1]*0.4*0.999;
    cdata=[linspace(1,F_color(1),Nm+1)',linspace(1,F_color(2),Nm+1)',linspace(1,F_color(3),Nm+1)'];
    cdata=reshape(cdata,Nm+1,1,3);
    if k>3
        patch([CX3_B(n:k),NaN],[CY3_B(n:k),NaN],1:Nm+1,'EdgeColor','interp','Marker','none',...
          'MarkerFaceColor','flat','CData',cdata,'LineWidth',1);
    end
    pause(0.02)
    %这里可以插入gif输出程序
end

function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn    ,yn       ,Input);
    K2=Fdydx(xn+h/2,yn+h/2*K1,Input);
    K3=Fdydx(xn+h/2,yn+h/2*K2,Input);
    K4=Fdydx(xn+h  ,yn+h*K3  ,Input);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
Output=[];
end

function dydx=Fdydx(x,y,Input)
%将原方程整理为dy/dx=F(y,x)的形式
%三摆方程
%输入Input整理
m=Input(2);
l=Input(3);
g=Input(4);
%输入
th=y(1:3);%角度1
pth=y(4:6);%角动量1
%利用拉格朗日法得到的方程
M=l^2*m*(-169+81*cos(2*(th(1)-th(2)))-9*cos(2*(th(1)-th(3)))+45*cos(2*(th(2)-th(3))));
dth1=(6*(-23*pth(1)+9*cos(2*(th(2)-th(3)))*pth(1)+27*cos(th(1)-th(2))*pth(2)-9*cos(th(1)+th(2)-2*th(3))*pth(2)+21*cos(th(1)-th(3))*pth(3)-27*cos(th(1)-2*th(2)+th(3))*pth(3)))/M;
dth2=(6*(27*cos(th(1)-th(2))*pth(1)-9*cos(th(1)+th(2)-2*th(3))*pth(1)-47*pth(2)+9*cos(2*(th(1)-th(3)))*pth(2)-27*cos(2*th(1)-th(2)-th(3))*pth(3)+57*cos(th(2)-th(3))*pth(3)))/M;
dth3=(6*(21*cos(th(1)-th(3))*pth(1)-27*cos(th(1)-2*th(2)+th(3))*pth(1)-27*cos(2*th(1)-th(2)-th(3))*pth(2)+57*cos(th(2)-th(3))*pth(2)-143*pth(3)+81*cos(2*(th(1)-th(2)))*pth(3)))/M;
dth=[dth1;dth2;dth3];
dpth1=-0.5*l*m*(5*g*sin(th(1))+l*dth(1)*(3*dth(2)*sin(th(1)-th(2))+dth(3)*sin(th(1)-th(3))));
dpth2=-0.5*l*m*(-3*l*dth(1)*dth(2)*sin(th(1)-th(2))+3*g*sin(th(2))+l*dth(2)*dth(3)*sin(th(2)-th(3)));
dpth3=0.5*l*m*(l*dth(1)*dth(3)*sin(th(1)-th(3))+l*dth(2)*dth(3)*sin(th(2)-th(3))-g*sin(th(3)));
%整理输出
dydx=zeros(6,1);
dydx(1)=dth1;
dydx(2)=dth2;
dydx(3)=dth3;
dydx(4)=dpth1;
dydx(5)=dpth2;
dydx(6)=dpth3;
end

3 多摆的混沌性

对于多摆来说，初始角度越小，其线性越好，非线性影响越小。比如下图为不同角度自由释放的双摆模型，随着时间推进，角度越大的摆，其与周边摆的差异性越大。在最下面的几个摆基本保持着左右摆动的周期运动。而靠上的摆早早的在第一个周期内，便出现了运动的混乱性（定性来看，就是率先的失去了和周围摆运动的一致性，代表误差经系统传递被放大到失去控制的程度）。

对于三摆来说，也是一样，越靠下线性度越好：

下面是三摆绘图的程序。双摆大同小异，稍微在基础上改一下就行，就不放出来了。

%三摆，做出很多摆的动画
%4阶RK求解
clear
clc
close all

%输入
N=3;%三摆
m=1;
l=1;
g=9.8;
Input=[N,m,l,g];
h=5e-3;
x0=0:h:20;
angle_list=1:2:359;
N_angle=numel(angle_list);
tN=numel(x0);
%构建储存单元
CX1_A_Sum=zeros(N_angle,tN);
CX1_B_Sum=CX1_A_Sum;
CX2_A_Sum=CX1_A_Sum;
CX2_B_Sum=CX1_A_Sum;
CX3_A_Sum=CX1_A_Sum;
CX3_B_Sum=CX1_A_Sum;

CY1_A_Sum=CX1_A_Sum;
CY1_B_Sum=CX1_A_Sum;
CY2_A_Sum=CX1_A_Sum;
CY2_B_Sum=CX1_A_Sum;
CY3_A_Sum=CX1_A_Sum;
CY3_B_Sum=CX1_A_Sum;
for k=1:N_angle
    angle_k=deg2rad(angle_list(k));
    %初始条件和时间设置
    y0=[angle_k;angle_k;angle_k;0;0;0];%这里全部是弧度值。分别代表[摆1与垂面夹角，摆2与垂面夹角，摆1角动量，摆2角动量]
    %代入到ODE求解器中
    [y1,Output]=ODE_RK4_hyh(x0,h,y0,Input);
    
    %提取出角度
    th1=y1(1,:);
    th2=y1(2,:);
    th3=y1(3,:);
    
    %计算出关节坐标
    CX1_A=zeros(1,tN);
    CX1_B=CX1_A+l*sin(th1);
    CY1_A=zeros(1,tN);
    CY1_B=CY1_A-l*cos(th1);
    
    CX2_A=CX1_B;
    CX2_B=CX2_A+l*sin(th2);
    CY2_A=CY1_B;
    CY2_B=CY2_A-l*cos(th2);
    
    CX3_A=CX2_B;
    CX3_B=CX3_A+l*sin(th3);
    CY3_A=CY2_B;
    CY3_B=CY3_A-l*cos(th3);
    
    CX1_A_Sum(k,:)=CX1_A;
    CX1_B_Sum(k,:)=CX1_B;
    CX2_A_Sum(k,:)=CX2_A;
    CX2_B_Sum(k,:)=CX2_B;
    CX3_A_Sum(k,:)=CX3_A;
    CX3_B_Sum(k,:)=CX3_B;
    
    CY1_A_Sum(k,:)=CY1_A;
    CY1_B_Sum(k,:)=CY1_B;
    CY2_A_Sum(k,:)=CY2_A;
    CY2_B_Sum(k,:)=CY2_B;
    CY3_A_Sum(k,:)=CY3_A;
    CY3_B_Sum(k,:)=CY3_B;
    
end

%绘图
%第t时刻的单摆图
t_k=100;
figure('Color','white')
set(gcf,'position',[488   342   400   300])
for t_k=1:4:800%720
clf
xlim([-N,N])
ylim([-N,N])
axis off equal
hold on
for k=1:N_angle
    %绘制摆
    plot([CX1_A_Sum(k,t_k),CX1_B_Sum(k,t_k)],[CY1_A_Sum(k,t_k),CY1_B_Sum(k,t_k)],'color','k')
    plot([CX2_A_Sum(k,t_k),CX2_B_Sum(k,t_k)],[CY2_A_Sum(k,t_k),CY2_B_Sum(k,t_k)],'color','k')
    plot([CX3_A_Sum(k,t_k),CX3_B_Sum(k,t_k)],[CY3_A_Sum(k,t_k),CY3_B_Sum(k,t_k)],'color','k')
end
hold off
pause(0.01)
end

function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn    ,yn       ,Input);
    K2=Fdydx(xn+h/2,yn+h/2*K1,Input);
    K3=Fdydx(xn+h/2,yn+h/2*K2,Input);
    K4=Fdydx(xn+h  ,yn+h*K3  ,Input);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
Output=[];
end

function dydx=Fdydx(x,y,Input)
%将原方程整理为dy/dx=F(y,x)的形式
%输入Input整理
m=Input(2);
l=Input(3);
g=Input(4);
%输入
th=y(1:3);%角度1
pth=y(4:6);%角动量1
%利用拉格朗日法得到的方程
M=l^2*m*(-169+81*cos(2*(th(1)-th(2)))-9*cos(2*(th(1)-th(3)))+45*cos(2*(th(2)-th(3))));
dth1=(6*(-23*pth(1)+9*cos(2*(th(2)-th(3)))*pth(1)+27*cos(th(1)-th(2))*pth(2)-9*cos(th(1)+th(2)-2*th(3))*pth(2)+21*cos(th(1)-th(3))*pth(3)-27*cos(th(1)-2*th(2)+th(3))*pth(3)))/M;
dth2=(6*(27*cos(th(1)-th(2))*pth(1)-9*cos(th(1)+th(2)-2*th(3))*pth(1)-47*pth(2)+9*cos(2*(th(1)-th(3)))*pth(2)-27*cos(2*th(1)-th(2)-th(3))*pth(3)+57*cos(th(2)-th(3))*pth(3)))/M;
dth3=(6*(21*cos(th(1)-th(3))*pth(1)-27*cos(th(1)-2*th(2)+th(3))*pth(1)-27*cos(2*th(1)-th(2)-th(3))*pth(2)+57*cos(th(2)-th(3))*pth(2)-143*pth(3)+81*cos(2*(th(1)-th(2)))*pth(3)))/M;
dth=[dth1;dth2;dth3];
dpth1=-0.5*l*m*(5*g*sin(th(1))+l*dth(1)*(3*dth(2)*sin(th(1)-th(2))+dth(3)*sin(th(1)-th(3))));
dpth2=-0.5*l*m*(-3*l*dth(1)*dth(2)*sin(th(1)-th(2))+3*g*sin(th(2))+l*dth(2)*dth(3)*sin(th(2)-th(3)));
dpth3=0.5*l*m*(l*dth(1)*dth(3)*sin(th(1)-th(3))+l*dth(2)*dth(3)*sin(th(2)-th(3))-g*sin(th(3)));
%整理输出
dydx=zeros(6,1);
dydx(1)=dth1;
dydx(2)=dth2;
dydx(3)=dth3;
dydx(4)=dpth1;
dydx(5)=dpth2;
dydx(6)=dpth3;
end