首页/文章/ 详情

求解滑动摆系统的运动形式

过冷水

3年前浏览3305

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，qq：927550334

一个单摆悬挂于可沿水平光滑轨道滑动的滑块上，如图8-1所示。滑块质量为m，单摆的摆杆长度为l，质量不计，摆锤质量为m'，摆杆与滑块光滑铰接，整个系统在同一竖直平面内运动，研究整个滑动摆系统的运动情况。

图 8-1 滑动摆模型示意图

如图8-1所示，沿直线水平轨道建立Ox轴，以x表示滑块在轨道上的位置，以θ表示摆杆与竖直线的夹角。则系统图 8-1 滑动摆模型示意图的拉格朗日函数为

求解式（8-2）可得滑动摆系统运动情况的解析表达式。

8.2 用 ode45 求解常微分方程

上述微分方程组比较复杂，用dsolve求解解析解非常耗时，且不一定可行。在实际应用中，对于复杂的微分方程（组），很多时候解析解是很难得到的，因此经常以求解数值解来代替解析解。本节介绍一个求解微分方程数值解的函数——ode45。ode45采用四阶和五阶Runge-Kutta单步算法，用变步长求解器求解非刚性常微分方程，其解具有二阶精度。在MATLAB中，ode45是解决微分方程（组）

句法

[T, Y] = ode45(odefun,tspan,y0,options)

说明

ode45：函数用来求解微分方程（组）数值的解。

odefun：函数句柄，可以是函数文件名、匿名函数句柄或内联函数名。

tspan：求解时间区间或时间点。

y0：初值向量。

options：求解参数（可选）,可以用 odeset 在计算前设定容许误差、输出参数、事件等。

T：返回列向量的时间点。

Y：返回对应 T 时刻的求解结果。

ode45求解常微分方程（组）的最重要一步就是对方程（组）进行降阶处理，使求解方程（组）全部变为一阶微分方程。对于式（8-2），可进行如下降阶处理。

进行降阶处理后，就可用ode45命令直接求解。为了更直观地了解滑动摆滑块与摆锤之间的相对运动，可以利用Maltab的可视化功能对滑动摆的运动进行可视化模拟。滑动摆的可视化过程与上一章傅科摆可视化类似，这里不再赘述，只给出程序。可视化程序完成两个功能：第一是绘制系统运动参数x和θ的时程变化曲线（图8-2），第二是生成一个模拟滑动摆运动过程的动画（avi视频文件），图8-3是其中一帧图像（动画制作的内容将在第9章介绍）。

图 8-2 滑块位置与摆线角度随时间变化曲线

图 8-3 某时刻滑动摆的运动状态图

%==================================================================% 
%% 文件名:SlidePendulum.m 
%% 功能:利用ode45求解滑动摆轨迹的数值解,将其可视化,同时将摆动过程制作成avi格式 
%% 的视频 
%==================================================================% 
% 
close all;clear all;clc 
g = 9.8;m1 = 3;m2 = 3;l = 1; 
[t, y] = ode45('SlidePendulum_fun',[0:0.05:5.5],... 
    [pi/4,0,-cos(pi/4)/2,0],[],m2/(m1 m2),g,l); 
%% 画曲线 
figure('color',[1 1 1],'unit','normalized','position',[0.1 0.1 0.45 
0.5]); 
[ax, h1, h2] = plotyy(t,y(:,3),t,180*y(:,1)/pi); 
xlabel('时间{\itt}','fontname','Times new roman','fontsize',20); 
ylabel(ax(1),'滑块位置{\itx}','fontname','Times new roman','fontsize', 
20)
ylabel(ax(2),'摆线角度{\it\theta}','fontname','Times new roman', 'fontsize', 
20) 
set(ax(1),'fontname','Times new roman','fontsize',20,'xlim',[0 5.5]) 
set(ax(2),'fontname','Times new roman','fontsize',20,'xlim',[0 5.5], 
'ylim',[-60 60]) 
set(h1,'linestyle','-','marker','*','linewidth',2) 
set(h2,'linestyle','-','marker','.') 
%% 可视化 
fig = figure('color',[1 1 1],'unit','normalized','position',[0.5 0.5 
0.45 0.45]); 
axis([-0.6 0.6 -1 0.2]);axis equal;axis off;hold on
y1 = -l*cos(y(:,1)); 
x1 = y(:,3) l.*sin(y(:,1)); 
line1 = line([-0.6,0.6],[0,0],'linewidth',4,'color',[0 0 0]);  
%滑轨初始化 
line2 = line([0,0],[0,-1],'linewidth',1,'linestyle','-.','color',[0 0 
0]); %竖轴初始化 
pole = line([y(1,3),x1(1)],[0,y1(1)],'color',[1 1 0],'linestyle', 
'-',...
'linewidth',3,'erasemode','xor'); %滑杆初始化 
block1 = line([y(1,3)-0.08,y(1,3) 0.08],[0,0],'color',[0 0 1],... 
    'linestyle','-','linewidth',20,'erasemode','xor');%滑块初始化 
block2 = line(x1(1),y1(1),'color',[1 0 0],'marker','.', 'markersize',... 
60,'erasemode','xor');%摆锤初始化 
mov = avifile('Video.avi','fps',20);  
% 创建视频文件 Video.avi,设帧率为 20fps 
for i = 1:size(t,1) 
    set(pole,'xdata',[y(i,3),x1(i)],'ydata',[0,y1(i)]); 
    set(block1,'xdata',[y(i,3)-0.08,y(i,3) 0.08],'ydata',[0,0]); 
    set(block2,'xdata',x1(i),'ydata',y1(i)); 
    drawnow; %更新滑动摆位置 
    pause(0.05) 
    F = getframe(fig); 
mov = addframe(mov,F); % 添加图像至视频文件 
end     
mov = close(mov); % 关闭视频文件
%==================================================================% 
%% 文件名 SlidePendulum_fun.m 
%% 本函数为上述程序中,计算灰度重心的子函数 
%% 参数 image:切割后的用于计算的图像 
%% 参数 threshold:灰度阈值  
%==================================================================% 
function k = SlidePendulum_fun(t,y,flag,M,g,l) 
k = [y(2); 
    (-M*sin(y(1))*cos(y(1))*y(2)^2-... 
    g/l*sin(y(1)))/(1-M*cos(y(1))^2); 
    y(4);
     (M*g*sin(y(1))*cos(y(1)) M*l*... 
    sin(y(1))*y(2)^2)/(1-M*cos(y(1))^2)];