Моделирование динамических систем: введение в GNU Octave

Большинство научных работников не ломают голову над тем, как устроены численные методы внутри. Они просто используют их, применяя в своей работе специализированные пакеты численных расчетов. Это совершенно не означает, что не нужно разбираться с тем, как эти методы устроены. Программу пишет человек, а ему свойственно ошибаться. И ошибки сквозят даже в самых дорогих и навороченных системах численной математики сплошь и рядом. К тому же есть задачи, где применение стандартных систем невозможно.

Вместе с тем, умение использовать универсальный математический софт это must have для современного ученого, потому что изобретая велосипед можно никогда не добраться до решения своей основной задачи. Сегодня мы рассмотрим обещанный Octave, попытавшись решить с его помощью очередную детскую задачу, сделав при этом недетские выводы.

1. Что такое Octave и где его взять

GNU Octave — свободный математический пакет численной математики, основанный на концепциях своего проприетарного собрата Matlab. Он доступен для всех популярных на сегодняшний день десктопных платформ, и получить его можно, зайдя в раздел загрузки на официальном сайте.

Большую часть экрана занимает так называемое командное окно, с приглашием командной строки вида ">>". Введем туда что нибудь и нажмем Enter

шикарно, переменная ans содержит результат последней операции, если вы явно не завели переменную под результат, например так

переменные вводятся «на ходу», как это принято в интерпретируемых языках. Да собственно m-язык октава и есть язык интерпретатора, очень напоминающий аналогичный язык матлаба. Один раз заведенная переменная сохраняет свое значение и тип до следующего присваивания, либо до полной очистки памяти командной оболочки

С левой стороны окна есть область, в которой расположены (сверху вниз): файловый менеджер, список глобальных переменных с их значениями, а также история введенных команд

Стоит обратить внимание на то, что для переменных указана размерность (Dimention) 1х1. Основным типом данных в Octave, как и в Matlab является матрица. Так вот наши числа это на самом деле матрицы размером 1 на 1 элемент.

Очистить все глобальные переменные можно выбрав в меню Edit -> Clear Workspace, а очистить командное окно — щелкнуть по нему правой кнопкой, выбрав Clear Window. Выполним очистку переменных и введя x получим

видно, что переменная x стала недоступна, до тех пока мы явно не присвоили ей значение снова.

Скажу честно, что с Ocatve я никогда не работал и разбираюсь вместе с читателем. Поэтому, предположим, что мы уже знаем достаточно, чтобы решить какую нибудь задачу. Например ту, которую решали в прошлый раз. Но сначала…

2. Нормальная форма Коши

и получаем вместо одного уравнения 2-го порядка два уравнения 1-го порядка, которые нужно решать вместе, одновременно

Это уже не одно уравнение, а система дифференциальных уравнений, содержащая две неизвестные зависимости и .

Как решать систему уравнений численно? Да так же, как и одно уравнение. Соберем все наши переменные в вектор-столбец (или массив, или матрицу-столбец, кому как удобно выражаться это не важно)

и правые части тоже соберем в вектор-столбец

Тогда рекуррентная формула метода Эйлера будет справедлива для каждого элемента этих столбцов

где n — число уравнений, в нашем случае их два.

К чему я всё это? А к тому, что Octave требует представлять дифференциальные уравнения именно в нормальной форме Коши.

С точки зрения механики вектор-столбец y называется вектором состояния системы или вектором фазовых координат. Тогда наша система уравнений превращается в одно векторное уравнение

Теперь возьмем, и в командном окне Octave наберем

Таким образом мы определяем функцию, которая на вход принимает текущее значение фазовых координат y и текущее время t, а на выходе возвращает значение производных от фазовых координат. Видим, что приглашение командной строки ">>" пропало, система ждет что мы закончим ввод тела функции ключевым словом endfunction. Продолжим набирать функцию

Итак, в теле функции мы определили g = 10 — принятое нами значение ускорения свободного падения. Мы видим, что же не появилось в списке переменных слева — эта переменная является локальной, и существует лишь в пределах функции. Переменная y это матрица-столбец, первый элемент y(1) которой это координата z, второй y(2) — проекция вектора скорости vz. Соответственно, dydt — это то значение, которое возвращает наша функция, это тоже матрица-столбец первый элемент которой это производная от z по времени, а второй — производная от vz по времени. То есть мы записали нашу систему дифференциальных уравнений в терминах Octave.

Теперь определимся с диапазоном времени, для которого нам надо получить результат. Пусть это будет десять точек от 0 до 1 секунды, с шагом 0.1 — сравним результат с ручным примером

Тут всё очевидно: t0 — начальный момент времени, tend — конечный момент времени. А вот deltaT = 0.1 секунда — это не шаг интегрирования! Это интервал, с которым Ocatve будет отображать для нас численное решение уравнения. Последняя команда формирует массив моментов времени для которых мы хотим получить решение.

Как вы знаете из прошлой статьи, для решения уравнения численно необходимо знать начальные значения скорости и координаты. Необходимо задать их для Octave

тем самым мы определили матрицу-столбец, содержащую начальную координату и начальное значение вертикальной проекции скорости. А теперь решаем уравнение

Функция lsode решит для нас уравнение численно. Первый параметр — имя функции, которая вычисляет производные от фазовых координат. Это функция F, мы её задали. Второй параметр — начальные условия, то есть значения фазовых координат в момент времени t = t0. Последний параметр — массив моментов времени, для которых мы хотим рассчитать значения фазовых координат. Жмем Enter…

В ответ вываливается гора требухи, оканчивающаяся предложением

нам предлагают листать требуху дальше (f), вернуться назад (b), или выйти (q). Те, кто знает *nix-подобные системы, знают, что это консольный вывод под управлением юниксовой утилиты less. Мы притворимся что не знаем, выйдем отсюда, нажав на клавиатуре q.

Наберем теперь в командной строке «y» и нажмем Enter

Ничего не напоминает? Ну конечно же это то самое решение, что мы получили для задачи из прошлой статьи. Только оно теперь удивительно точное — значения совпадают с аналитическим решением! Открою вам тайну — это точное решение нашей задачи. Связано это с тем, что функция lsode использует для решения задачи отнюдь не метод Эйлера, а нечто более продвинутое. Движение камня происходит с постоянным ускорением, и та формула численной аппроксимации, что применяется в данном методе, очевидно просто совпадает с аналитическим решением задачи. Хотя, если лезть в дебри представления машиной чисел с плавающей запятой, то… А ну да ладно, сейчас не об этом.

3. Строим график

Наберите теперь команду

Выскочит окошко с графическим представлением решения

Синяя кривая сверху — координата камня. Оранжевая снизу — вертикальная проекция скорости. Введенная команда неудобна тем, что строит график для всех функций сразу. А если мы хотим, построить, зависимость скорости от высоты? Тогда поступим так

Будет построен график, где по оси абсцисс пойдет переменная y(1), а по оси ординат — y(2), что есть соответственно высота и вертикальная проекция скорости

Такой график в механике называют фазовым портретом системы, то есть траекторией движения системы в пространстве фазовых координат. В данном случае фазовые координаты это высота камня над поверхностью Земли и его скорость. В механике часто нельзя получить точного решения в виде функций времени, но можно связать между собой координаты и скорости исходя из энергетических соображений. По полученным зависимостям можно найти фазовую траекторию системы и сделать много выводов о характере её движения не решая задачу до конца. Об этом мы тоже как-нибудь поговорим.

А пока предлагаю выполнить самостоятельное задание — постройте график зависимости высоты от времени и скорости от времени в отдельных окнах. И старайтесь не смотреть под спойлер

Заключение

На этом первое знакомство с Octave, как средой для численного моделирования будем считать оконченным. Теперь у нас есть мощный инструмент, с помощью которого мы будем постигать более изощренные премудрости.