Main Content

이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.

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

이 예제에서는 공중으로 던져진 바통의 동특성을 모델링하는 연립상미분방정식을 풉니다 [1]. 바통은 길이가 L인 막대로 연결되고 질량이 m1m2인 두 개의 입자로 모델링됩니다. 바통은 공중으로 던져지고 이후 중력으로 인한 힘에 따라 수직 xy 평면에서 이동합니다. 막대는 수평과의 각도 θ를 형성하고 첫 번째 질량의 좌표는 (x,y)입니다. 이 정식화에서 두 번째 질량의 좌표는 (x+L cos θ,y+L sin θ)입니다.

세 개의 좌표 x, y, θ에 각각 라그랑주 방정식을 적용하면 연립 운동 방정식이 얻어집니다.

(m1+m2)x¨-m2Lθ¨sinθ-m2Lθ˙2cosθ=0,(m1+m2)y¨-m2Lθ¨cosθ-m2Lθ˙2sinθ+(m1+m2)g=0,L2θ¨-Lx¨sinθ+Ly¨cosθ+gLcosθ=0.

MATLAB®에서 이 연립 ODE를 풀려면 솔버 ode45를 호출하기 전에 방정식을 함수로 코딩하십시오. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.

방정식 코딩하기

ode45 솔버를 사용하려면 방정식을 q˙=f(t,q) 형식으로 작성해야 합니다. 여기서 q˙는 각 좌표의 1계 도함수입니다. 이 문제에서 해 벡터에는 6개의 구성요소가 있습니다. x, y, 각도 θ 및 각각의 1계 도함수입니다.

q=[q1q2q3q4q5q6]=[xx˙yy˙θθ˙].

이 표기법을 사용하여 운동 방정식을 q의 요소로만 다시 작성할 수 있습니다.

(m1+m2)q2˙-m2Lq6˙sinq5=m2Lq62cosq5,(m1+m2)q4˙-m2Lq6˙cosq5=m2Lq62sinq5-(m1+m2)g,L2q6˙-Lq2˙sinq5+Lq4˙cosq5=-gLcosq5.

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

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

[1000000m1+m2000-m2Lsinq5001000000m1+m20m2Lcosq50000100-Lsinq50Lcosq50L2][q1˙q2˙q3˙q4˙q5˙q6˙]=[q2m2Lq62cosq5q4m2Lq62sinq5-(m1+m2)gq6-gLcosq5]

이 표현식에서 질량 행렬의 0이 아닌 요소를 계산하는 함수를 작성할 수 있습니다. 함수는 세 개의 입력을 받습니다. 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) q˙=f(t,q)의 각 방정식 우변에 대한 함수를 작성할 수 있습니다. 질량 행렬 함수와 마찬가지로 이 함수는 tq에 대한 두 개의 필수 입력과 파라미터 값 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

연립방정식 풀기

먼저 구조체에서 적절하게 명명된 필드를 설정하여 m1, m2, 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);

문제의 초기 조건을 설정합니다. 바통은 일정 각도로 위로 던져졌으므로 초기 속도에 0이 아닌 값을 사용합니다(x0˙=4, y0˙=20, θ0˙=2). 초기 위치의 경우 바통은 수직으로 세워진 위치에서 시작하므로 θ0=-π/2, x0=0, y0=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. The axes with title Motion of a Thrown Baton, Solved by ODE45 contains 75 objects of type line.

참고 문헌

[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.

로컬 함수(Local Function)

여기 나열된 함수는 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

참고 항목

관련 항목