Mathematica 模拟滑块摆运动

写这个主要是为了Mathematica的动画制作，Mathematica官方文档给出的是用Animate及相关的函数，但这并不能满足题目的要求，所以选择创建GIF图片。

## 题目

滑块摆由一置于光滑杆上的质量为m的滑块A、一质量为M的小球B和长度为L，质量不计的刚性杆铰接而成，不计各处摩擦，以过A点的水平面为零势能面，通过Lagrange方程建立系统的运动方程，利用Maple软件画出：
1. 滑块A的位移x随时间t的变化曲线
2. 角度φ随时间t的变化曲线
3. 滑块摆的运动动画

![](http://blog.qiniu.ctaoist.cn/滑块摆示意图.png)

## 拉格朗日力学分析
在力学中经常使用 Lagrange 方程来分析。
在滑块摆这个简单的物理模型中，假如用牛顿第二定律，则必须仔细地辨明所有的相关作用力。这是一项既困难又容易出错的工作。
如下定义的拉格朗日方程:

$$
\begin{align}
L &= E_c-E_p \\
\frac{d}{dt}\left(\frac{\partial L}{\partial q'}\right)-\frac{\partial L}{\partial q}+\frac{\partial D}{\partial q'} &= F_q
\end{align}
$$

其中:
 - q x(t)和θ(t)的自由度
 - D 由于摩擦而消耗的能量
 - Fq 由自由度 q 产生的力
 - Ec 和 Ep 系统的动能和势能

系统有两个自由度，以 x 和 $\varphi$ 为广义坐标，以过 A 点的水平面为零势能面，
系统的动能和势能分别为:

$$
\begin{align}
E_c &= \frac{1}{2}m\dot{x}^2+\frac{1}{2}M(\dot{x}^2+l^2\dot{\varphi}^2+2l\dot{x}\dot{\varphi}\cos{\varphi}) \newline
&= \frac{1}{2}(m+M)\dot{x}^2+\frac{1}{2}Ml^2\dot{\varphi}^2+Ml\dot{x}\dot{\varphi}\cos{\varphi} \newline
E_p &= -Mgl\cos{\varphi}
\end{align}
$$

所以，系统的 Lagrange 方程为

$$
\begin{align}
L&=E_c-E_p
=\frac{1}{2}(m+M)\dot{x}^2+\frac{1}{2}Ml^2\dot{\varphi}^2+Ml\dot{x}\dot{\varphi}\cos{\varphi}+Mgl\cos{\varphi}
\end{align}
$$

计算出诸导数

$$
\begin{align}
\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right) &= (m+M)\ddot{x}+Ml\ddot{\varphi}\cos{\varphi}-Ml\dot{\varphi}^2\sin{\varphi} \\
\frac{\partial L}{\partial x} &= 0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\varphi}}\right) &= Ml^2\ddot{\varphi}+Ml\ddot{x}\cos{\varphi}-Ml\dot{\varphi}\dot{x}\sin{\varphi} \\
\frac{\partial L}{\partial \varphi} &= -Ml\dot{\varphi}\dot{x}\sin{\varphi}-Mgl\sin{\varphi}
\end{align}
$$

带入 Lagrange 方程，得到系统的运动微分方程

$$
\begin{align}
(m+M)\ddot{x}+Ml\ddot{\varphi}\cos{\varphi}-Ml\dot{\varphi}^2\sin{\varphi}=0 \\
l\ddot{\varphi}+\ddot{x}\cos{\varphi}+g\sin{\varphi}=0
\end{align}
$$

## Mathematica源码

*某些代码可在[Mma脚本](http://mmascript.com/)在线执行，这份代码不行，因为含有`Export`等被禁用的命令。*
```mathematica
Clear[x, \[Phi], fig];
m = M = 1; L = 2; g = 9.8;
solve = NDSolve[{
    (m + M) x''[t] + M L \[Phi]''[t] Cos[\[Phi][t]] - 
      M L \[Phi]'[t]^2 Sin[\[Phi][t]] == 0,
    L \[Phi]''[t] + x''[t] Cos[\[Phi][t]] + g Sin[\[Phi][t]] == 0,
    \[Phi][0] == 0, x[0] == 0, \[Phi]'[0] == -1.3, x'[0] == 1},
   {x, \[Phi]},
   {t, 0, 20}
   ];
{x, \[Phi]} = {x, \[Phi]} /. solve[[1]];
Plot[x[t], {t, 0, 20}, AxesLabel -> {t, "x"}]
Plot[\[Phi][t], {t, 0, 20}, AxesLabel -> {t, "\[Phi]"}]
fig = {};
Do[
 huakuai = {x[t], 0}; 
 ball = {x[t] + L Sin[\[Phi][t]], -L Cos[\[Phi][t]]};(*滑块和球的坐标*)
 
 tmpFig = Show[Graphics[Disk[huakuai, 0.1], PlotLabel -> {黑色滑块, 红色为球}],
   Graphics[{Red, Disk[ball, 0.1]}],
   Graphics[Line[{huakuai, ball}]], Axes -> {True, False}, 
   PlotRange -> {{-7, 1}, {-3, 1}}];
 AppendTo[fig, tmpFig], {t, 0, 20, 0.1}]
Export["i:/data/滑块摆的运动动画.gif", fig](*将图画列表导出为gif*)
```

## 动画结果

![](http://blog.qiniu.ctaoist.cn/滑块摆的运动动画.gif)

Mathematica 模拟滑块摆运动

最近更新

热门标签

文章目录