1. Модель
2. Функция правых частей
3. Запуск процесса численного интегрирования
4. Передача в функцию правых частей дополнительных параметров
5. Управления точностью численного интегрирования

Модель

Рассмотрим простейшее дифференциальное уравнение математического маятника \[\ddot{\varphi} = - \frac{g}{l} \sin \varphi\]

где \(g\) - ускорение свободного падения, \(l\) - длина подвеса материальной точки, \(\varphi\) – угол отклонения маятника от вертикали.

Необходимо численно проинтегрировать это уравнение для следующих начальных условий: \[\varphi(0) = 60^{\circ}, \quad \dot{\varphi}(0) = 0.\]

Для того, чтобы в MATLAB численно проинтегрировать дифференциальное уравнения второго порядка, это уравнение необходимо привести к системе двух дифференциальных уравнений первого порядка. Для этого введём новую переменную – угловую скорость \(\omega = \dot{\varphi}\), и запишем следующую эквивалентную систему дифференциальных уравнений: \[\begin{cases} \dot{\varphi} = \omega, \\ \dot{\omega} = - \frac{g}{l} \sin \varphi. \end{cases}\]

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

Полученную систему дифференциальных уравнений можно представить в матричной форме: \[\dot{\mathbf{q}} = \mathbf{f}(t, \mathbf{q})\]

где \(\mathbf{q}\) – вектор состояния системы (матрица-столбец): \[\mathbf{q} = \begin{bmatrix} \varphi \\ \omega \end{bmatrix}\]

\(\mathbf{f}\) – векторная функция векторного аргумента и, может быть, времени. Эту функцию далее будем называть функцию правых частей системы дифференциальных уравнений: \[\mathbf{f}(t, \mathbf{q}) = \begin{bmatrix} \omega \\ - \frac{g}{l} \sin \varphi \end{bmatrix} = \begin{bmatrix} q_2 \\ - \frac{g}{l} \sin q_1 \end{bmatrix}\]

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

function dq = f(t, q)
    % Для наглядности дадим элементам вектора состояния
    % более понятные имена,  
    phi = q(1); %  1 элемент q - угол
    w   = q(2); %  2 элемент q - угловая скорость

    g   = 9.807; % ускорение свободного падения
    l   = 1.0;   % длина подвеса груза

    % Первое уравнение системы диф. уравнений
    % Производная фи = угловая скорость
    dphi  = w;
    % Второе уравнение системы
    % Производная угловой скорости = угловое ускорение
    dw    = -g*sin(phi)/l;

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

Запуск процесса численного интегрирования

Для запуска процесса численного интегрирования MATLAB предлагает несколько функций, отличающихся методом используемого численного интегрирования (ode45, ode113, ode23, …). Например, для интегрирования методом Рунге-Кутты 4-го и 5-го порядков с автоматическим выбором шага интегрирования в зависимости от заданной точности используется функция ode45.

Создадим и запустим файл-скрипт main.m со следующим кодом:

% Вектор-столбец начальных условия [угол; угловая скорость]
q0 = [60*pi/180; 0];

% Интервал интегрирования от 0 до 2 с
ts = [0, 2];

% Запускаем процесс интегрирования:
% Вызываем функцию ode45, передавая ей ссылку на функцию правых частей
% интервал интегрирования и вектор-столбец начальных условий
[t, q] = ode45(@f, ts, q0);

Вектор начальных условий q0 содержит начальный угол поворота маятника (первый элемент) и угловую скорость (второй элемент). Этот же порядок следования переменных состояния ожидается в функции правых частей.

Функция ode45 в этом примере записывает результат своей работы в столбец t и прямоугольную матрицу q. В столбце t содержатся моменты времени, для которых получены решения. Первый элемент в этом столбце равен 0 с, последний 2 c. Количество промежуточных значений алгоритм численного интегрирования определяет автоматически, исходя из требований точности интегрирования.

Матрица q содержит результаты интегрирования. Число строк этой матрицы равно числу строк столбца t, а в каждой строке матрицы q(k,:) содержится вектор-строка состояния системы, соответствующий моменту времени t(k).

Далее в файл-скрипте main.m можно построить график зависимости от времени угла отклонения маятника:

plot(t, q(:,1)*180/pi);
xlabel('t, c');
ylabel('\phi, градус');

Передача в функцию правых частей дополнительных параметров

В написанной выше функции правых частей параметры рассматриваемой механической системы \(g\) и \(l\) были заданы непосредственно в коде функции, что нельзя назвать хорошей практикой программирования, поскольку смешивается код и данные модели. Отделим код модели от её данных, для этого перенесем определение всех параметров системы в главный файл-скрипт main.m, а сами параметры будем передавать в функцию правых частей третьим аргументом. Новое определение функции будет иметь следующий вид:

function dq = f(t, q, p)
    % Для наглядности дадим элементам вектора состояния 
    % более понятные имена
    phi = q(1); %  1 элемент вектора q - угол
    w   = q(2); %  2 элемент вектора q - угловая скорость

    % Для сокращения записи уравнений движения
    % перепишем значения полей структуры 
    % в новые локальные переменные
    g = p.g;
    l = p.l;

    % Первое уравнение системы
    dphi  = w;
    % Второе уравнение системы
    dw    = -g*sin(phi)/l;

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

В функцию правых частей третьим аргументом передаётся структура p, с атрибутами (полями) p.l и p.g, значения которых используются в коде функции. Структуру необходимо создать в главном файл-скрипте main.m:

% Параметры системы
p.g = 9.807; % ускорение свободного падения
p.l = 1.0;   % длина подвеса

% Начальные условия
q0 = [60*pi/180; 0];
% Интервал интегрирования
ts = [0, 2];

Функция ode45 ожидает, что первым аргументов будет ссылкой (@) на функцию двух переменных (время и вектор состояния), но новая функция правых частей это функция трех переменных – f(t,q,p). Для адаптации новой функции f используем анонимную функцию-обёртку двух переменных (t, q), которая будет вызывать функцию f, передавая в нее t и q, а также структуру с параметрами системы.

f2 = @(t, q) f(t,q,p);
[t, q] = ode45(@f2, ts, q0);

или в одну строчку

[t, q] = ode45(@(t, q) f(t,q,p), ts, q0);

Управления точностью численного интегрирования

Для управления точностью интегрирования в функцию ode45 необходимо передать дополнительный параметр (структуру), которая создается при помощи специальной функции odeset. Например, для того чтобы задать относительную погрешность равную \(1 \cdot 10^{-8}\), а абсолютную погрешность \(1 \cdot 10^{-10}\) необходимо вызвать функцию со следующими аргументами:

options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);

При численном интегрировании MATLAB будет пытаться выполнить условия \[\varepsilon_i \leq \max(\text{RelTol} \cdot \max(\text{abs}(q)), \text{AbsTol}),\]

для всех переменных состояния, самостоятельно изменяя шаг интегрирования, где \(\varepsilon_i\) – относительна погрешность вычисления \(i\)-ой переменной состояния.

Точность можно задать для каждой переменной отдельно:

options = odeset('RelTol',[1e-7, 1e-9],'AbsTol', [1e-9, 1e-8]);

Здесь мы задали относительную погрешность \(10^{-7}\) для угла и \(10^{-9}\) для угловой скорости, и абсолютную погрешность \(10^{-9}\) для угла и \(10^{-8}\) для угловой скорости. В этом случае для каждой \(i\)-ой переменной состояния будет выполнятся условие: \[\varepsilon_i \leq \max(\text{RelTol}_i \cdot \text{abs}(q_i)), \text{AbsTol}_i)\]

Для запуска численного интегрирования с заданной точностью необходимо передать структуру options в функцию ode45 последним аргументом:

[t, q] = ode45(@(t, q) f(t,q,p), ts, q0, options);