Решите уравнения движения для маркера, выданного в воздух

Скрипт Open 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.

Уравнения кода

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 как структура, которая содержит значения параметров и затем использует метод, описанный в Параметризации Функций, чтобы передать структуру функции как дополнительный вход.

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] Уэллс, отважьтесь схему А. Шаума теории и проблемы лагранжевой динамики: с обработкой уравнений Эйлером движения, уравнений Гамильтона и принципа Гамильтона. Ряд схемы Шаума. Нью-Йорк: паб Schaum. Ко, 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

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