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

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

С помощью этой записи можно переписать уравнения движения полностью в терминах элементов $\mathit{q}$:

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

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

Из этого выражения можно записать функцию, которая вычисляет ненулевые элементы массовой матрицы. Функция принимает три входа: $\mathit{t}$ и $\mathit{}$требуется вектор решения q (необходимо указать эти входы, даже если они не используются в теле функции), а $\mathit{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

Далее можно записать функцию для правой части каждого из уравнений в системе $\mathit{M}\mathit{(}t\mathit{,}\text{q}\stackrel{}{\mathit{)}}\mathit{}\mathit{q\dot{}=f}\mathit{(}t$, q). Как и функция матрицы масс, эта функция принимает два обязательных входа $\mathit{}$для t $\mathit{}$и q и один необязательный вход для $\mathit{}$значений параметров 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 значений параметров для ${\mathit{}}_{\mathrm{m1}}$, ${\mathit{}}_{\mathrm{m2}}$, $\mathit{g}$и $\mathit{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);

Задайте начальные условия проблемы. Так как дубинка брошена вверх под углом, используйте ненулевые значения для начальных скоростей, $\stackrel{}{{\mathit{}}_{\mathrm{}}}x0\dot{}=4\mathrm{,}$$\stackrel{}{{\mathit{}}_{\mathrm{}}}\mathrm{y0\dot{}=20}$и $\stackrel{}{{}_{\mathrm{}}}\theta 0\dot{}=2\mathrm{.}$Для исходных положений полицейская дубинка начинается в вертикальном положении, таким образом${\mathrm{,\; \theta 0}}_{\mathrm{}}=-\pi /2\mathrm{,}$${\mathit{}}_{\mathrm{}}x0=0\mathrm{,}$и ${\mathit{}}_{\mathrm{}}y0=L\mathit{.}$

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

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

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

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

[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

Ссылки

[1] Уэллс, Dare A. Schaum наброски теории и проблемы лагранжевой динамики: с лечением уравнений движения Эйлера, уравнений Гамильтона и Гамильтона. Серия набросков Шаума. Нью-Йорк: Шаум Паб. Ко, 1967.

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

Здесь перечислены локальные вспомогательные функции, вызываемые решателем 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