Моделирование динамических систем: введение в GNU Octave16.02.2018 15:03

Жили-были умные, но очень жадные люди, которые написали замечательную программу Matlab. Умные они были потому, что программа вышла хорошей, а жадными, потому что очень любили деньги. Так любили, что брали их за свой Matlab не только с дядек серьезных, матлабом деньги зарабатывающих, а и с бедных студентов тоже, которым порой и сухую корочку хлеба купить не за что было. И кончилась бы сказочка скоро и невесело, если бы мир был не без добрых и умных людей, написавших похожие на матлаб программы, хоть худо-бедно работающие, да для всех желающих бесплатные. И с открытыми исходными текстами. Так что сами бедные студенты стали те программы дописывать, и работать они лучше и лучше стали с каждым годом. И стали тогда все жить-поживать, да добра наживать…

Введение

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

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

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

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

>> 2 + 2
ans = 4

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

>> x = 3 + 2
x =  5

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

>> x + ans
ans =  9

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

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

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

>> x
error: 'x' undefined near line 1 column 1
>> x = 3
x =  3
>> x
x =  3

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

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

2. Нормальная форма Коши
Помните, я говорил о том, что большинство, за исключением некоторого узкого класса, численных методов «заточены» под уравнения 1-го порядка? И что прежде чем применить эти методы к решению уравнения, то надо понизить его порядок? Беда в том, что не все уравнения допускают понижение порядка. Однако, любое уравнение n-го порядка можно превратить в n уравнений 1-го порядка. Это называется преобразованием уравнения к нормальной форме Коши. Возьмем прошлое уравнение

$\ddot z = -g$

Вспоминаем, что

$\dot z = v_z$

тогда

$\dot v_z = -g$

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

$\begin{cases} &\dot z = &v_z \\ &\dot v_z = &-g \end{cases}$

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

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

$\mathbf y = \begin{bmatrix} z \\ v_z \end{bmatrix}$

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

$\mathbf F(t, \mathbf y) = \begin{bmatrix} v_z \\ -g \end{bmatrix}$

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

$y_j^{(i+1)} = y_j^{(i)} + F_j^{(i)} \, \Delta t, \quad j = \overline{1,n}$

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

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

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

$\frac{d\mathbf y}{dt} = \mathbf F(\mathbf y, t)$

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

>> function dydt = F(y, t)

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

>> function dydt = F(y, t)

g = 10

dydt(1) = y(2)
dydt(2) = -g

endfunction

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

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

>> t0 = 0
t0 = 0
>> tend = 1
tend =  1
>> deltaT = 0.1
deltaT =  0.10000
>> t = [t0:deltaT:tend]
t =

   0.00000   0.10000   0.20000   0.30000   0.40000   0.50000   0.60000   0.70000   0.80000   0.90000   1.00000

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

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

>> y0 = [100; 0]
y0 =

   100
     0

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

>> y = lsode("F", y0, t)

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

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

g =  10
dydt = -0.99690
dydt =

   -0.99690  -10.00000

g =  10
dydt = -0.99690
dydt =

   -0.99690  -10.00000

g =  10
dydt = -0.99690
dydt =

   -0.99690  -10.00000

g =  10
dydt = -1.9936
-- less -- (f)orward, (b)ack, (q)uit

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

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

>> y
y =

   100.00000     0.00000
    99.95000    -1.00000
    99.80000    -2.00000
    99.55000    -3.00000
    99.20000    -4.00000
    98.75000    -5.00000
    98.20000    -6.00000
    97.55000    -7.00000
    96.80000    -8.00000
    95.95000    -9.00000
    95.00000   -10.00000

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

3. Строим график
Наберите теперь команду

>> plot(t, y)

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

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

>> plot(y(:,1), y(:,2))

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

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

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

Ответ

Высота z (t)

>> plot(t, y(:,1))

Скорость vz (t)

>> plot(t, y(:,2))

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

Спасибо за внимание, увидимся!