求解抛向空中的短棒的运动方程
此示例求解常微分方程组,该方程组对抛向空中的短棒的动态进行建模 [1]。短棒建模为两个质点 和 由长度为 的棒连接。短棒抛向空中,随后在重力作用下在垂直 xy 平面内移动。棒与水平线构成角度 ,第一个质点的坐标是 。根据此公式,第二个质点的坐标是 。
系统的运动方程是通过对三个坐标 、 和 中的每个坐标应用拉格朗日方程获得的:
编写方程代码
MATLAB 要求将方程写作 形式,其中 是每个坐标的一阶导数。在此问题中,解向量有六个分量:、、角度 及各自的一阶导数:
通过此表示法,您可以用 的元素彻底重写运动方程:
遗憾的是,运动方程不符合求解器要求的形式 ,因为左侧有几个具有一阶导数的项。出现这种情况时,必须使用质量矩阵来表示方程的左侧。
通过矩阵表示法,您可以使用 形式的质量矩阵将运动方程重写为包含六个方程的方程组。质量矩阵用矩阵-向量积表示方程左侧的一阶导数的线性组合。在此形式中,方程组变为:
根据此表达式,您可以编写一个用于计算质量矩阵的非零元素的函数。该函数接受三个输入: 和解向量 (这两个输入是必需的,即使这些输入在函数体中不使用,您也必须指定这些输入),以及 (可选的额外输入,用于传入参数值)。
function M = mass(t,q,P) % Extract parameters m1 = P(1); m2 = P(2); L = P(3); g = P(4); % Mass matrix elements 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
接下来,您可以为方程组 中的每个方程的右侧编写一个函数。与质量矩阵函数一样,此函数接受两个必需的输入 和 ,以及一个用于参数值的可选输入 。
function dydt = f(t,q,P) % Extract parameters m1 = P(1); m2 = P(2); L = P(3); g = P(4); % Equations 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
求解方程组
首先,为 、、 和 创建一个由参数值组成的向量 P
。求解器将向量 P
作为额外输入传递给质量矩阵和 ODE 函数。
P = [0.1 0.1 1 9.81];
为问题创建一个初始条件向量。由于短棒以一定角度向上投掷,因此初始速度应使用非零值,即 、 和 。短棒的初始位置从直立位置开始,即 、 和 。
y0 = [0; 4; P(3); 20; -pi/2; 2];
现在,创建一个 ode
对象来表示问题,为 ODEFcn
、MassMatrix
、InitialValue
和 Parameters
指定值。对象显示表明,因为存在质量矩阵,已为该问题选择 ode15s
求解器。
F = ode(ODEFcn=@f,... MassMatrix=@mass,... InitialValue=y0,... Parameters=P)
F = ode with properties: Problem definition ODEFcn: @f Parameters: [0.1000 0.1000 1 9.8100] InitialTime: 0 InitialValue: [6x1 double] Jacobian: [] MassMatrix: [1x1 odeMassMatrix] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: ode15s Show all properties
为积分时间跨度创建一个由 0 和 4 之间的 25 个点组成的向量,并使用 solve
仿真方程组。当您指定时间向量时,求解器会采用它自己的内部步长,但在您指定的点对解进行插值。
tspan = linspace(0,4,25); sol = solve(F,tspan)
sol = ODEResults with properties: Time: [0 0.1667 0.3333 0.5000 0.6667 0.8333 1 1.1667 1.3333 1.5000 1.6667 1.8333 2 2.1667 2.3333 2.5000 2.6667 2.8333 3 3.1667 3.3333 3.5000 3.6667 3.8333 4] Solution: [6x25 double]
绘制结果
solve
的输出包含运动方程在每个请求的时间步的解。为了检验结果,对短棒随时间的运动轨迹绘图。
遍历解的每列,并在每个时间步绘制短棒的位置。给短棒的每端涂上不同颜色,以便看到它随时间的旋转。
figure title("Motion of Thrown Baton"); axis([0 22 0 25]) hold on for j = 1:length(sol.Time) theta = sol.Solution(5,j); X = sol.Solution(1,j); Y = sol.Solution(3,j); xvals = [X X+P(3)*cos(theta)]; yvals = [Y Y+P(3)*sin(theta)]; plot(xvals,yvals,xvals(1),yvals(1),"ro",xvals(2),yvals(2),"go") end hold off
参考
[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.