求解抛向空中的短棒的运动方程

打开实时脚本

此示例求解常微分方程组，该方程组对抛向空中的短棒的动态进行建模 [1]。短棒建模为两个质点 $m_{1}$ 和 $m_{2}$ 由长度为 $L$ 的棒连接。短棒抛向空中，随后在重力作用下在垂直 xy 平面内移动。棒与水平线构成角度 $θ$ ，第一个质点的坐标是 $(x, y)$ 。根据此公式，第二个质点的坐标是 $(x + L \cos θ, y + L \sin θ)$ 。

系统的运动方程是通过对三个坐标 $x$ 、 $y$ 和 $θ$ 中的每个坐标应用拉格朗日方程获得的：

$\begin{array}{l} (m_{1} + m_{2}) \ddot{x} - m_{2} L \ddot{θ} \sin θ - m_{2} L {\dot{θ}}^{2} \cos θ = 0, \\ (m_{1} + m_{2}) \ddot{y} - m_{2} L \ddot{θ} \cos θ - m_{2} L {\dot{θ}}^{2} \sin θ + (m_{1} + m_{2}) g = 0, \\ L^{2} \ddot{θ} - L \ddot{x} \sin θ + L \ddot{y} \cos θ + g L \cos θ = 0 . \end{array}$

要在 MATLAB® 中求解此 ODE 方程组，请在调用求解器 ode45 之前将这些方程编码为一个函数。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

ode45 求解器要求将方程写作 $\dot{q} = f (t, q)$ 形式，其中 $\dot{q}$ 是每个坐标的一阶导数。在此问题中，解向量有六个分量： $x$ 、 $y$ 、角度 $θ$ 及各自的一阶导数：

$q = [\begin{array}{c} q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ q_{5} \\ q_{6} \end{array}] = [\begin{array}{c} x \\ \dot{x} \\ y \\ \dot{y} \\ θ \\ \dot{θ} \end{array}] .$

通过此表示法，您可以用 $q$ 的元素彻底重写运动方程：

$\begin{array}{l} (m_{1} + m_{2}) \dot{q_{2}} - m_{2} L \dot{q_{6}} \sin q_{5} = m_{2} L q_{6}^{2} \cos q_{5}, \\ (m_{1} + m_{2}) \dot{q_{4}} - m_{2} L \dot{q_{6}} \cos q_{5} = m_{2} L q_{6}^{2} \sin q_{5} - (m_{1} + m_{2}) g, \\ L^{2} \dot{q_{6}} - L \dot{q_{2}} \sin q_{5} + L \dot{q_{4}} \cos q_{5} = - g L \cos q_{5} . \end{array}$

遗憾的是，运动方程不符合求解器要求的形式 $\dot{q} = f (t, q)$ ，因为左侧有几个具有一阶导数的项。出现这种情况时，必须使用质量矩阵来表示方程的左侧。

通过矩阵表示法，您可以使用 $M (t, q) \dot{q} = f (t, q)$ 形式的质量矩阵将运动方程重写为包含六个方程的方程组。质量矩阵用矩阵-向量积表示方程左侧的一阶导数的线性组合。在此形式中，方程组变为：

$[\begin{array}{c} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & m_{1} + m_{2} & 0 & 0 & 0 & - m_{2} L \sin q_{5} \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & m_{1} + m_{2} & 0 & m_{2} L \cos q_{5} \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & - L \sin q_{5} & 0 & L \cos q_{5} & 0 & L^{2} \end{array}] [\begin{array}{c} \dot{q_{1}} \\ \dot{q_{2}} \\ \dot{q_{3}} \\ \dot{q_{4}} \\ \dot{q_{5}} \\ \dot{q_{6}} \end{array}] = [\begin{array}{c} q_{2} \\ m_{2} L q_{6}^{2} \cos q_{5} \\ q_{4} \\ m_{2} L q_{6}^{2} \sin q_{5} - (m_{1} + m_{2}) g \\ q_{6} \\ - g L \cos q_{5} \end{array}]$

根据此表达式，您可以编写一个用于计算质量矩阵的非零元素的函数。该函数接受三个输入： $t$ 和解向量 $q$ （这两个输入是必需的，即使这些输入在函数体中不使用，您也必须指定这些输入），以及 $P$ （可选的额外输入，用于传入参数值）。要将此问题的参数传递给函数，请创建 P 作为保留参数值的结构体，然后使用参数化函数中所述的方法将该结构体作为额外输入传递给函数。

function M = mass(t,q,P)
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Mass matrix function
M = zeros(6,6);
M(1,1) = 1;
M(2,2) = m1 + m2;
M(2,6) = -m2*L*sin(q(5));
M(3,3) = 1;
M(4,4) = m1 + m2;
M(4,6) = m2*L*cos(q(5));
M(5,5) = 1;
M(6,2) = -L*sin(q(5));
M(6,4) = L*cos(q(5));
M(6,6) = L^2;
end

接下来，您可以为方程组 $M (t, q) \dot{q} = f (t, q)$ 中的每个方程的右侧编写一个函数。与质量矩阵函数一样，此函数接受两个必需的输入 $t$ 和 $q$ ，以及一个用于参数值的可选输入 $P$ 。

function dydt = f(t,q,P)
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Equation to solve
dydt = [q(2)
        m2*L*q(6)^2*cos(q(5))
        q(4)
        m2*L*q(6)^2*sin(q(5))-(m1+m2)*g
        q(6)
        -g*L*cos(q(5))];
end

求解方程组

首先，通过在结构体中设置适当命名的字段，为 $m_{1}$ 、 $m_{2}$ 、 $g$ 和 $L$ 创建参数值结构体 P。结构体 P 作为额外的输入传递给质量矩阵和 ODE 函数。

P.m1 = 0.1;
P.m2 = 0.1;
P.L = 1;
P.g = 9.81

P = struct with fields:
    m1: 0.1000
    m2: 0.1000
     L: 1
     g: 9.8100

为积分时间跨度创建一个由 0 和 4 之间的 25 个点组成的向量。当您指定时间向量时，ode45 会在请求的时间返回插值解。

tspan = linspace(0,4,25);

设置问题的初始条件。由于短棒以一定角度向上投掷，因此初始速度应使用非零值，即 $\dot{x_{0}} = 4$ 、 $\dot{y_{0}} = 20$ 和 $\dot{θ_{0}} = 2$ 。短棒的初始位置从直立位置开始，即 $θ_{0} = - π / 2$ 、 $x_{0} = 0$ 和 $y_{0} = L$ 。

y0 = [0; 4; P.L; 20; -pi/2; 2];

使用 odeset 创建一个引用质量矩阵函数的 options 结构体。由于求解器预期质量矩阵函数只有两个输入，使用匿名函数从工作区传入参数 P。

opts = odeset('Mass',@(t,q) mass(t,q,P));

最后，使用 ode45 求解方程组，输入如下：

方程 @(t,q) f(t,q,P) 的匿名函数
请求解的时间组成的向量 tspan
初始条件组成的向量 y0
引用质量矩阵的 options 结构体 opts

[t,q] = ode45(@(t,q) f(t,q,P),tspan,y0,opts);

绘制结果

ode45 的输出包含运动方程在每个请求的时间步的解。为了检验结果，对短棒随时间的运动轨迹绘图。

遍历解的每行，并在每个时间步绘制短棒的位置。给短棒的每端涂上不同颜色，以便看到它随时间的旋转。

figure
title('Motion of a Thrown Baton, Solved by ODE45');
axis([0 22 0 25])
hold on
for j = 1:length(t)
   theta = q(j,5);
   X = q(j,1);
   Y = q(j,3);
   xvals = [X X+P.L*cos(theta)];
   yvals = [Y Y+P.L*sin(theta)];
   plot(xvals,yvals,xvals(1),yvals(1),'ro',xvals(2),yvals(2),'go')
end
hold off

Figure contains an axes object. The axes object with title Motion of a Thrown Baton, Solved by ODE45 contains 75 objects of type line. One or more of the lines displays its values using only markers

参考

[1] Wells, Dare A. Schaum’s Outline of Theory and Problems of Lagrangian Dynamics:With a Treatment of Euler’s Equations of Motion, Hamilton’s Equations and Hamilton’s Principle.Schaum's Outline Series.New York:Schaum Pub.Co, 1967.

局部函数

此处列出 ODE 求解器为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

function dydt = f(t,q,P) % Equations to solve
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% RHS of system of equations
dydt = [q(2)
        m2*L*q(6)^2*cos(q(5))
        q(4)
        m2*L*q(6)^2*sin(q(5))-(m1+m2)*g
        q(6)
        -g*L*cos(q(5))];
end
%------------------------------------------------
function M = mass(t,q,P) % Mass matrix function
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Set nonzero elements in mass matrix
M = zeros(6,6);
M(1,1) = 1;
M(2,2) = m1 + m2;
M(2,6) = -m2*L*sin(q(5));
M(3,3) = 1;
M(4,4) = m1 + m2;
M(4,6) = m2*L*cos(q(5));
M(5,5) = 1;
M(6,2) = -L*sin(q(5));
M(6,4) = L*cos(q(5));
M(6,6) = L^2;
end

另请参阅

ode45