求解滑动摆系统的运动形式
过冷水诚挚邀请你加入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是解决微分方程（组）
上述微分方程组比较复杂，用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)];
过冷水发表于仿真秀 平台原创文章，未经授权禁止私自转载，如需转载请需要和作者沟通表明授权声明，未授权文章皆视为侵权行为，必将追责。如果您希望加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载均可加群：927550334。
精品回顾
过冷水和你分享 matlab读取存储各种文件的方法 文末有独家金曲分享
