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

  • Tutorial
Жили-были умные, но очень жадные люди, которые написали замечательную программу 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, как средой для численного моделирования будем считать оконченным. Теперь у нас есть мощный инструмент, с помощью которого мы будем постигать более изощренные премудрости.

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

Комментарии 16

    0
    Продукт хороший, но после matlab смотрится не очень. Ну и и симулинк это замечательное решение, аналогов которому нет.

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

      Ничего подобного
      ПК МВТУ
      SciLab
        0

        Я не знаю, что там сделали МВТУ, хотя рискну предположить, что это вряд ли является полноценным аналогом для Simulink.
        Но вот SciLab является аналогом только на уровне несложных (студенческих) работ. Как чуть что посложнее, и посыпалось. У меня аспиарнт был, сторонник open-source, он хотел всю кандидатскую сделать только с использованием Octave и SciLab. И если рассчёты в Octave ещё туда-сюда, то от моделирования в SciLab пришлось отказаться, не взлетело.

          0
          Ну и что, что не является.
          Помнится и линукс ну совсем не являлся аналогом Windows NT на серверах. И что? Какова доля этого ядра на серверах и в топ500 суперкомпьютеров?

          Пусть не является. Пусть. Зато он есть
            0
            Пусть не является, пусть падает на собственных примерах. Главное что open-source.
              –1
              Главное что open-source

              Да!
                +1
                А то что падает, так перестанет падать со временем. И обсвятится это без тех, кто постоянно по этому поводу ноет
                  0
                  Вот когда починят, посмотрим. А пока 55 usd для студента лучшее решение, имхо, если ему конечно это вообще нужно.
              0
              Ну порой самописный код получше готовых пакетов будет, это сильно зависит от задачи.
            0
            А чем matlab лучше? Подробности.
              +1

              Вот здесь неплохо сказано:
              The MathWorks has many toolboxes for MATLAB, there's Simulink and its related products for which there really is no equivalent in Octave (yes, you'd have to pay for all that. But often your employer/school does that anyway, and well, it at least exists), proven compliance with several industry standards, testing tools, validation tools, requirement management systems, report generation, a much larger community & user base, etc. etc. etc. MATLAB is only a small part of something much larger. Octave is...just Octave.

                +1
                Вот знаете что, по-моему я нигде не сказал, что Octave полностью заменяет матлаб. И вообще, статья не об этом. Вам не нравится Octave? Не используйте Octave. Вам не нравится SciLab? Не используйте SciLab. Вам Linux не нравится? Не используйте Linux.

                Только вот, в конце 90-х, когда я начал заниматься программированием, чтобы начать нужно было:

                1. Купить компьютер.
                2. Зайти на пиратский развал и купить там диск с ОС и интересующей средой и инструментарием. В сумме всё это понянет не на одну сотню тысяч рублей.

                Сегодня нужно выполнить лишь первый пункт, купить железо. Всё остальное можно взять даром. Вам от этого плохо? Плохо от того что не нужно кланяться охреневшим в конец от жадности Билли, Стиву и прочим? Вам нужен мир в котором, для того чтобы спеть песню военных лет на празднике Победы нужно заплатить господину Михалкову? Если вас это устраивает — ради бога.

                Только вот оффтоп попрошу прекратить.
                0
                Octave тормозит. В смысле считает медленнее даже 2010 Матлаба (там, где время существенно)
              0

              После примера уровня x=3+2 уже можно отойти от командной строки в сторону создания m-файла.
              ПК МВТУ — интересная штука, но рассматривать ее как аналог Симулинка все-таки слишком, хоть и для определенного круга задач они сходны.

                +1
                Будет m-файл, куда торопится. Статьи намеренно небольшого объема для того чтобы мозги не запутывать
                0
                Octave все же неполноценная замена MATLAB: местами синтаксис разный, а в версии 3.8 (в Debian) нет LDL- декомпозиции, например (хотя, может, я не так понял help).
                На фоне python+scipy и R я не вижу особого смысла в использовании Octave.

                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                Самое читаемое