煤化工网

MATLAB非线性可视化(引3)多摆模型

网友投稿
后台-系统设置-扩展变量-手机广告位-内容正文顶部

MATLAB非线性可视化(引3)多摆模型

原创 hyhhyh21 matlab爱好者
matlab爱好者

matlabaihaozhe

学matlab编程就关注matlab爱好者!

收录于话题

接着前面的Mandelbrot集和牛顿迭代继续介绍一个非线性模型:多摆。如果只看到前面的两个引子,肯定会有疑问:非线性只是一种通过迭代产生的数学游戏吗?
事实上,非线性存在于物理与工程中的各个领域。在机械中,就存在着大量的非线性现象。通过双摆和三摆的例子,来感受到一个小的扰动,随着时间的推移,到最终会带来多大的变化。
这里不涉及到计算多体动力学和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上发布

标签:

后台-系统设置-扩展变量-手机广告位-内容正文底部
留言与评论(共有 0 条评论)
   
验证码: