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

Этот пример решает систему обыкновенных дифференциальных уравнений, которые моделируют динамику маркера, выданного в воздух [1]. Маркер моделируется как две частицы с массами m1 и m2 соединенный стержнем длины L. Маркер выдан в воздух и впоследствии перемещается в вертикальную xy-плоскость, удовлетворяющую силе из-за силы тяжести. Стержень формирует угол θ с горизонталью и координатами первой массы (x,y). С этой формулировкой координаты второй массы (x+Lcosθ,y+Lsinθ).

Уравнения движения для системы получены путем применения уравнений Лагранжа для каждой из трех координат, 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®, закодируйте уравнения в функцию прежде, чем вызвать решатель ode45. Можно или включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

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

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

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.

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

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

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

От этого выражения можно записать функцию, которая вычисляет ненулевые элементы большой матрицы. Функция берет три входных параметров: 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). Как функция большой матрицы, эта функция берет два необходимых входных параметров для 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 из значений параметров для m1, m2, 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);

Установите начальные условия проблемы. Поскольку маркер выдан вверх под углом, используйте ненулевые значения для начальных скоростей, x0˙=4, y0˙=20, и θ0˙=2. Для исходных положений маркер начинается в вертикальном положении, таким образом, θ0=-π/2, x0=0, и y0=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

Ссылки

[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

Смотрите также

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте