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

라이브 스크립트 열기

이 예제에서는 공중으로 던져진 바통의 동특성을 모델링하는 연립상미분방정식을 풉니다[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]
         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

Figure contains an axes object. The axes object with title Motion of Thrown Baton 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.

참고 항목

ode45