공중으로 던져진 바통의 운동 방정식 풀기
이 예제에서는 공중으로 던져진 바통의 동특성을 모델링하는 연립상미분방정식을 풉니다[1]. 바통은 길이가 인 막대로 연결되고 질량이 과 인 두 개의 입자로 모델링됩니다. 바통은 공중으로 던져지고 이후 중력으로 인한 힘에 따라 수직 xy 평면에서 이동합니다. 막대는 수평과의 각도 를 형성하고 첫 번째 질량의 좌표는 입니다. 이 정식화에서 두 번째 질량의 좌표는 입니다.
세 개의 좌표 , , 에 각각 라그랑주 방정식을 적용하면 연립 운동 방정식이 얻어집니다.
방정식 코딩하기
MATLAB에서는 방정식을 형식으로 작성해야 합니다. 여기서 는 각 좌표의 1계 도함수입니다. 이 문제에서 해 벡터에는 6개의 성분이 있습니다. , , 각도 및 각각의 1계 도함수입니다.
이 표기법을 사용하여 운동 방정식을 의 요소로만 다시 작성할 수 있습니다.
그러나 이 운동 방정식은 좌변에 1계 도함수가 있는 여러 항이 있기 때문에 솔버에서 요구하는 형식 에 맞지 않습니다. 이 경우 질량 행렬을 사용하여 방정식의 좌변을 표현해야 합니다.
행렬 표기법을 사용하여 질량 행렬을 사용한 6개의 연립방정식으로 운동 방정식을 형식으로 다시 작성할 수 있습니다. 질량 행렬은 행렬-벡터 곱을 사용하여 방정식의 좌변에 1계 도함수의 선형 결합을 표현합니다. 이 형식에서 연립방정식은 다음과 같이 됩니다.
이 표현식에서 질량 행렬의 0이 아닌 요소를 계산하는 함수를 작성할 수 있습니다. 함수는 세 개의 입력을 받습니다. 와 해 벡터 는 필수 입력이며(함수 본문에 사용되지 않더라도 이 입력을 지정해야 함) 는 파라미터 값으로 전달하는 데 사용되는 선택적 추가 입력입니다.
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];
문제의 초기 조건으로 구성된 벡터를 만듭니다. 바통은 일정 각도로 위로 던져졌으므로 초기 속도에 0이 아닌 값을 사용합니다(, , ). 초기 위치의 경우 바통은 수직으로 세워진 위치에서 시작하므로 , , 입니다.
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
참고 문헌
[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.