MATLAB非线性可视化（引3）多摆模型
事实上，非线性存在于物理与工程中的各个领域。在机械中，就存在着大量的非线性现象。通过双摆和三摆的例子，来感受到一个小的扰动，随着时间的推移，到最终会带来多大的变化。
这里不涉及到计算多体动力学和SimMechanics，只采用简单的拉格朗日方法建立运动模型。详细可见参考资料[1]。
双摆方程可以用两个摆的角度进行描述，其运动的角速度可以用角动量来描述。这样，我们就有了双摆的4个方程：
这个方程组为一个4阶的常微分方程，利用经典的龙格库塔RK方法，就可以进行求解。
随着时间的推移，双摆的运动越来越无序，变得难以琢磨。
如果把很多双摆在不同角度同时释放，可以得到如下的图像：
可以看到越靠上的双摆运动越混乱，靠近下方的双摆运动几乎像单摆一样。
对于三摆来说，方程也可以通过拉格朗日方法建立，但是结果会复杂一些。同样，我们还是定义变量为三个角度和三个角动量。这里方程比较繁琐，就直接放在程序里了。
最终的效果图如下：
同样，将多个三摆在不同角度同时释放的结果如下图：
当然，根据方程还可以求得更多阶的摆，然而计算量会越来越大。所以对于更高阶的摆，就不适用于直接求解出方程的方法了。
下面是双摆的代码。用到了自己编写的龙格库塔方法。
%双摆
%4阶RK求解
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:1160
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)
F=getframe(gca);
I=frame2im(F);
[I,map]=rgb2ind(I,20);
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再下面是三摆的程序。方程和双摆类似，但是由于增加了一个摆导致代码更长。
%三摆
%4阶RK求解
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)
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参考资料：
[1]https://en.wikipedia.org/wiki/Double_pendulum
图片来源：由 在Pixabay上发布
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2021-10-13，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读