Решите уравнения движения для Батона, брошенного в воздух

Открыть Live Script

Этот пример решает систему обыкновенных дифференциальных уравнений, которые моделируют динамику палочки, выброшенной в воздух [1]. Дубинка моделируется как две частицы с массами $m_{1}$ и $m_{2}$ соединенный стержнем длины $L$ . Дубинка выбрасывается в воздух и впоследствии перемещается в вертикальной 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 ®, кодируйте уравнения в функцию перед вызовом решателя ode45. Можно либо включить необходимые функции в качестве локальных функций в конец файла (как это сделано здесь), либо сохранить их как отдельные именованные файлы в директории по пути MATLAB.

Кодовые уравнения

The ode45 решатель требует, чтобы уравнения были записаны в виде $\dot{q} = f (t, q)$ , где $\dot{q}$ является первой производной от каждой координаты. В этой задаче вектор решения имеет шесть компонентов: $x$ , $y$ , угол $θ$ и их первые производные:

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

К сожалению, уравнения движения не вписываются в форму $\dot{q} = f (t, q)$ требуется решателем, так как слева несколько членов с первыми производными. Когда это происходит, необходимо использовать большую матрицу, чтобы представлять левую сторону уравнения.

С помощью матричного обозначения можно переписать уравнения движения как систему шести уравнений с помощью большой матрицы в виде $M (t, q) \dot{q} = f (t, q)$ . Большая матрица выражает линейные комбинации первых производных в левой части уравнения с матрицей-векторного продукта. В этом виде система уравнений становится:

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

Из этого выражения можно написать функцию, которая вычисляет ненулевые элементы большой матрицы. Функция принимает три входов: $t$ и вектор решения $q$ требуются (необходимо задать эти входы, даже если они не используются в теле функции), и $P$ является необязательным дополнительным входом, используемым для передачи значений параметров. Чтобы передать параметры для этой задачи функции, создайте P как структура, которая содержит значения параметров, а затем использует метод, описанный в Parameterizing Functions, чтобы передать структуру в функцию в качестве дополнительного входа.

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) \dot{q} = f (t, q)$ . Как и функция большой матрицы, эта функция принимает два необходимых входов для $t$ и $q$ и один необязательный вход для значений параметров $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

Решающая система уравнений

Во-первых, создайте структуру P значений параметров для $m_{1}$ , $m_{2}$ , $g$ , и $L$ путем установки соответствующих именованных полей в структуре. Структура P передается в большую матрицу, и ОДУ функционирует как дополнительный вход.

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

Создайте вектор с 25 точками от 0 до 4 для временного интервала интегрирования. Когда вы задаете вектор времени, ode45 возвращает интерполированные решения в требуемое время.

tspan = linspace(0,4,25);

Установите начальные условия задачи. Поскольку дубинка выброшена вверх под углом, используйте ненулевые значения для начальных скоростей, $\dot{x_{0}} = 4$ , $\dot{y_{0}} = 20$ , и $\dot{θ_{0}} = 2$ . Для начальных позиций дубинка начинается в вертикальном положении, поэтому $θ_{0} = - π / 2$ , $x_{0} = 0$ , и $y_{0} = L$ .

y0 = [0; 4; P.L; 20; -pi/2; 2];

Использование odeset чтобы создать структуру опций, которая ссылается на функцию большой матрицы. Поскольку решатель ожидает, что функция большой матрицы будет иметь только два входов, используйте анонимную функцию, чтобы пройти в параметрах P из рабочей области.

opts = odeset('Mass',@(t,q) mass(t,q,P));

Наконец, решите систему уравнений, используя ode45 с этими входами:

Анонимная функция для уравнений @(t,q) f(t,q,P)
Векторная tspan времени, когда запрашивается решение
Векторная y0 начальных условий
Структура опций 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: С обработкой уравнений движения Эйлера, уравнений Гамильтона и принципа Гамильтона. Серия контуров Шаума. Нью-Йорк: Шаум Паб. Ко, 1967.

Локальные функции

Здесь перечислены локальные вспомогательные функции, которые вызывается решатель ОДУ для вычисления решения. Также можно сохранить эти функции как собственные файлы в директории по пути 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

См. также

ode45

Документация