공중으로 던져진 바통의 운동 방정식 풀기

라이브 스크립트 열기

이 예제에서는 공중으로 던져진 바통의 동특성을 모델링하는 연립상미분방정식을 풉니다[1]. 바통은 길이가 $L$ 인 막대로 연결되고 질량이 $m_{1}$ 과 $m_{2}$ 인 두 개의 입자로 모델링됩니다. 바통은 공중으로 던져지고 이후 중력으로 인한 힘에 따라 수직 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에서는 방정식을 $\dot{q} = f (t, q)$ 형식으로 작성해야 합니다. 여기서 $\dot{q}$ 는 각 좌표의 1계 도함수입니다. 이 문제에서 해 벡터에는 6개의 성분이 있습니다. $x$ , $y$ , 각도 $θ$ 및 각각의 1계 도함수입니다.

$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}$

그러나 이 운동 방정식은 좌변에 1계 도함수가 있는 여러 항이 있기 때문에 솔버에서 요구하는 형식 $\dot{q} = f (t, q)$ 에 맞지 않습니다. 이 경우 질량 행렬을 사용하여 방정식의 좌변을 표현해야 합니다.

행렬 표기법을 사용하여 질량 행렬을 사용한 6개의 연립방정식으로 운동 방정식을 $M (t, q) \dot{q} = f (t, q)$ 형식으로 다시 작성할 수 있습니다. 질량 행렬은 행렬-벡터 곱을 사용하여 방정식의 좌변에 1계 도함수의 선형 결합을 표현합니다. 이 형식에서 연립방정식은 다음과 같이 됩니다.

$[\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}]$

이 표현식에서 질량 행렬의 0이 아닌 요소를 계산하는 함수를 작성할 수 있습니다. 함수는 세 개의 입력을 받습니다. $t$ 와 해 벡터 $q$ 는 필수 입력이며(함수 본문에 사용되지 않더라도 이 입력을 지정해야 함) $P$ 는 파라미터 값으로 전달하는 데 사용되는 선택적 추가 입력입니다.

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

다음으로 연립방정식 $M (t, q) \dot{q} = f (t, q)$ 의 각 방정식 우변에 대한 함수를 작성할 수 있습니다. 질량 행렬 함수와 마찬가지로 이 함수는 $t$ 및 $q$ 에 대한 두 개의 필수 입력과 파라미터 값 $P$ 에 대한 하나의 선택적 입력을 받습니다.

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

연립방정식 풀기

먼저 $m_{1}$ , $m_{2}$ , $L$ , $g$ 에 대한 파라미터 값으로 구성된 벡터 P를 만듭니다. 솔버는 벡터 P를 질량 행렬 및 ODE 함수에 추가 입력으로 전달됩니다.

P = [0.1 0.1 1 9.81];

문제의 초기 조건으로 구성된 벡터를 만듭니다. 바통은 일정 각도로 위로 던져졌으므로 초기 속도에 0이 아닌 값을 사용합니다( $\dot{x_{0}} = 4$ , $\dot{y_{0}} = 20$ , $\dot{θ_{0}} = 2$ ). 초기 위치의 경우 바통은 수직으로 세워진 위치에서 시작하므로 $θ_{0} = - π / 2$ , $x_{0} = 0$ , $y_{0} = L$ 입니다.

y0 = [0; 4; P(3); 20; -pi/2; 2];

이번에는 ODEFcn, MassMatrix, InitialValue, Parameters의 값을 지정하여 문제를 나타내는 ode 객체를 만듭니다. 객체 표시를 보면 질량 행렬이 있기 때문에 문제에 대해 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]

   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.

참고 항목

ode45