Интегрирование уравнений движения эллиптического мятника в MATLAB

Рассмотрим задачу численного интегрирования уравнений движения эллиптического маятника, схема которого показана на рисунке 1. Эллиптический маятник состоит из поступательно движущейся платформы 1 и груза 2, подвешенного на невесомом стержне, который шарнирно закрепляется на платформе в точке А. Известны масса платформы \(m_1\), масса груза \(m_2\) и длина подвеса \(l = AB\).

Рисунок 1 - Эллиптический маятник

Положение системы определяется перемещением платформы относительно некоторой точки O на горизонтальной прямой (обобщенная координата \(x\)), и углом отклонения подвеса от вертикальной прямой (\(\varphi\)). Движение системы происходит в поле силы тяжести. Ускорение свободного падения направлено вдоль оси \(y_0\).

Уравнения движения эллиптического маятника имеют вид: \[\left\{ \begin{aligned} & \ddot x \cos \varphi + \ddot \varphi l = -g \sin \varphi \\ & \ddot x (m_1+m_2) + \ddot \varphi l m_2 \cos \varphi = l m_2 \dot{\varphi}^2 \sin \varphi \end{aligned} \right.\]

Для численного интегрирования этой системы дифференциальных уравнений в MATLAB необходимо:

• преобразовать её к системе дифференциальных уравнений первого порядка, введя дополнительные переменные: \(v = \dot x\) - линейная скорость платформы; \(\omega = \dot \varphi\) - угловая скорость поворота подвеса,
• привести систему к форме Коши, выразив производные вектора состояния системы \((x, \varphi, v, \omega)\):

Правые части последних двух уравнений системы могут быть получены из исходной системы двух дифференциальных уравнений второго порядка, которые представляют собой систему линейных уравнения относительно вторых производных: \[\left\{ \begin{aligned} & a_{11} \ddot x + a_{12} \ddot \varphi = b_1 \\ & a_{21} \ddot x + a_{22} \ddot \varphi = b_2 \end{aligned} \right.\]

или в матричной форме \[\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} \ddot x \\ \ddot \varphi \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}\]

где \[\begin{aligned} & a_{11} & = & \cos \varphi \\ & a_{12} & = & l \\ & a_{21} & = & m_1 + m_2 \\ & a_{22} & = & l m_2 \cos \varphi \\ & b_{1} & = & - g \sin \varphi \\ & b_{2} & = & l m_2 \dot{\varphi}^2 \sin \varphi \end{aligned}\]

В результате решения этой системы линейных уравнений находятся значения \(\ddot x, \ddot \varphi\) \[\ddot x = \dot v = f_v(t, x, \varphi, v, \omega)\]

и \[\ddot \varphi = \dot \omega = f_\omega(t, x, \varphi, v, \omega)\]

В рассматриваемом случае необходимо решить систему двух линейных уравнений, что можно сделать аналитически и полученные выражения \(f_v(t, x, \varphi, v, \omega)\) и \(f_\omega(t, x, \varphi, v, \omega)\) использовать при разработке файл-функции правых частей дифференциальных уравнений.

Для систем с большим числом степеней свободы аналитическое решение системы линейных относительно вторых производных обобщенных координат может быть очень громоздким, поэтому целесообразно решать эту систему численно в файл-функции правых частей.

Файл-функция правых частей

Создадим файл-функцию правых частей dq_elliptic(t, q, p), которая принимает в качестве аргументов текущее время, вектор состояния системы (столбец) и структуру с её параметрами. Примем следующий порядок обобщенных координат и скоростей в векторе q: \(x\), \(\varphi\), \(v\), \(\omega\).

function dq = dq_elliptic(t, q, p)
    % Даем имена элементам q для удобства чтения кода 
    x   = q(1);
    phi = q(2);
    v   = q(3);
    w   = q(4);

    % Формируем матрицу коэффициентов aij (построчно)
    A = [cos(phi), p.l; p.m1+p.m2, p.m2*p.l*cos(phi)];
    % Формируем матрицу-столбец B
    B = [-p.g*sin(phi); p.l*p.m2*w*w*sin(phi)];

    % Решаем систему линейных уравнений 
    % относительно вторых производных 
    % обобщенных координат
    accel = A\B;

    % Извлекаем из решения
    % ускорение платформы 
    ax = accel(1);
    % и угловое ускорение поворота стержня
    dw = accel(2);

    % формируем результат работы функции -- производную вектора состояния
    dq = [v; w; ax; dw];

Главный файл-скрипт

% Параметры системы
p.l  = 1;
p.m1 = 5;
p.m2 = 2;
p.g  = 9.807;

% Начальные условия
x0   = 0;
phi0 = 60*pi/180;
v0   = 0;
w0   = 0; 

q0 = [x0; phi0; v0; w0];

% Интервал интегрирования
% Получить решение на интервале от 0 до 5 секунд с шагом 0.01 с 
ts = 0:0.01:5;

% Запускаем процесс интегрирования
[t, q] = ode45(@(t,q) dq_elliptic(t, q, p), ts, q0);

% Графики перемещения платформы и угла поворота от времени 
subplot(2,1, 1);
plot(t, q(:,1));
xlabel('t, c');
ylabel('x, м');

subplot(2,1, 2);
plot(t, q(:,2));
xlabel('t, c');
ylabel('\phi, градус');

Матрица масс

Задачу разрешения уравнений относительно старших производных (решение системы линейных уравнений относительно производных) можно делегировать MATLAB-у исключив её из кода файл-функции правых частей. Для этого необходимо создать отдельную файл-функцию, которая вычисляет (возвращает) матрицу коэффициентов при производных системы дифференциальных уравнений.

Для рассматриваемого примера систему дифференциальных уравнений первого порядка в матричной форме можно записать в виде: \[\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & a_{11} & a_{12} \\ 0 & 0 & a_{21} & a_{22} \\ \end{bmatrix} \begin{bmatrix} \dot x \\ \dot \varphi \\ \dot v \\ \dot \omega \\ \end{bmatrix} = \begin{bmatrix} v \\ \omega \\ b_1 \\ b_2 \\ \end{bmatrix}\]

Матрица коэффициентов при производных имеет размерность \(4 \times 4\). Создадим файл-функцию mass.m, которая вычисляет эту матрицу для заданного момента времени и вектора состояния системы. Третьим аргументом также будем передавать в эту функцию структуру с параметрами системы:

function m = mass(t, q, p)
    phi = q(1);    
    m = [1,  0,  0,   0;
         0,  1,  0,   0;
         0,  0,  p.l, cos(phi); 
         0,  0,  p.l*p.m2*cos(phi), (p.m1+p.m2)];    
end

Файл-функция правых частей (right_side.m) будет возвращать только матрицу-столбец правой части системы: \[\begin{bmatrix} v \\ \omega \\ b_1 \\ b_2 \\ \end{bmatrix}\]

% Функция правых частей
function dq = right_side(t,q,p)
    phi = q(1);
    w   = q(3);            
    dq = [q(3); q(4); -p.g*sin(phi); p.l*p.m2*w*w*sin(phi)];
end

В главном файл-скрипте перед запуском процесса интегрирования необходимо сформировать структуру с дополнительными параметрами для интегратора при помощи функции odeset.

Атрибуту ‘Mass’ необходимо присвоить ссылку на файл-функцию, вычисляющую “матрицу масс” (это должна быть функция двух переменных, поэтому используем анонимную функцию-“обертку”).

При вызове функции ode45 первым аргументом будет ссылка на функцию, вычисляющую правую часть системы дифференциальных уравнений – right_side. Последним аргументом передается структура opt с параметрами интегрирования:

...
...

% Сообщаем численному методу дополнительные параметры, в 
% том числе ссылку на функцию, вычисляющую матрицу коэффициентов 
% при производных mass(...)
opt = odeset('Mass', @(t,q) mass(t,q,p), 'RelTol', 1e-5);

% Запускаем процесс интегрирования
[t, q] = ode45(@(t,q) right_side(t, q, p), ts, q0, opt);

...
...