Решение систем обыкновенных дифференциальных уравнений в среде MATLAB. Выполнение скриптов в процессе решения и отслеживание событий. Часть 3

В предыдущей публикации мы рассматривали опции общего назначения, но не все. В стороне остались 2 опции (Events и OutputFcn). В этой публикации мы детально рассмотрим их синтаксис и особенности, посмотрим примеры и поразмышляем относительно их возможностей использования.

1. Опция OutputFcn позволяет выполнять какие-то команды или даже целые скрипты прямо в процессе расчета решения системы диффуров.
2. Опция Events позволяет отслеживать события.

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

1. Начнем с опции отслеживания событий. Для этого служит опция Events. Она позволяет отслеживать события по заданному условию. Давайте посмотрим, что же эта опция может делать. Сначала разберемся с синтаксисом и функционалом, а потом порассуждаем о применении событий. 

В качестве параметра этой опции передается функция, которая должна содержать всю информацию об отслеживаемом событии или о нескольких событиях. Обратите внимание, что используется освещённый нами в первой публикации по данной теме оператор “собака”.

Как построить функцию событий? На рисунке показан синтаксис функции событий с одним условием. 

Функция событий должна иметь как минимум 2 аргумента: 

1. текущее время или аргумент t,
2. текущий вектор компонент решения X. 

А вот возвращаемых значений должно быть 3: 

1. position, 
   isterminal,
   direction.

Рассмотрим их подробнее.

В position записывается выражение, которое равно нулю при совершении события. Например, если вы хотите отследить момент, когда 3-й компонент решения равен какому-нибудь числу C, то условие будет выглядеть так, как показано на рисунке. 

Второе выходное значение – isterminal – может принимать 2 значения: 0 или 1. Как можно догадаться из названия, если isterminal – 1, то после попадания на событие, расчет прекращается. Если isterminal – 0, то расчет, соответственно, продолжается. 

Третье выходное значение – direction – отвечает за динамику изменения значения position. Если direction – 1, то событие считается состоявшимся если position = 0 и само значение position при этом возрастало. Если direction = -1, то значение position должно убывать. Если direction = 0, то событие считается состоявшимся вне зависимости от возрастания или убывания position. Пояснительная картинка показана на рисунке, в первом случае direction = 1 и отслеживается только первое событие. Во втором, когда direction = -1, отслеживается только второе. И в третьем, при равенстве нулю – оба обнуления position считаются событиями. 

Мы задали события. Теперь хотелось бы получить значения решений в точках, соответствующих событиям?

В прошлой публикации мы узнали, как информация о решении системы дифференциальных уравнений записывается в виде структуры MATLAB. Эта структура имеет свои поля, о которых вы можете узнать из help-а.

Но нас в данном случае интересуют поля xe и ye, относящиеся к событиям:

• Xe – значение аргумента в момент события.
• Ye – вектор значений компонент в момент события.

Посмотрим как это выглядит в скрипте.

В качестве системы диффуров возьмём задачу движения шарика с сопротивлением воздуха. Задаём начальные условия, задаём интервал времени, решаем.

Добавим событие – обнуление координаты y, или “касание земли”. Задаём direction, расчет не останавливаем поэтому isterminal 0, и направление движения нам без разницы – поэтому direction 0. 

Отобразим на графике событие в виде точки.

Изменим isterminal на 1 и запустим скрипт.

Как видим, решение на этой точке обрывается. 

Мы рассмотрели отслеживание одного события. А как быть с несколькими?

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

Здесь 2 события: первое – обнуление вертикальной компоненты скорости – X(4), а второе событие – обнуление вертикальной координаты.

Два события отобразились. Теперь, например, остановим счет, когда случится второе событие, а когда случится первое – останавливать не будем. Изменим 2-ю компоненту isterminal с нуля на единицу. Запустим скрипт.

Работает. 

Далее посмотрим работу direction. Если поставим direction = 1 для второго события, то в нашем примере оно не будет отслеживаемым, так как значение вертикальной координаты снижается. Проверим это. Запустим скрипт.

Как видим, это событие пропущено. 

Вот так работает отслеживания событий. 

Остановка вычисления при выполнении определённых условий поможет Вам сократить машинное время, а также позволяет реализовывать что-то вроде примера, который вы сейчас увидите.

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

2. Выполнение скриптов в процессе вычисления

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

Эта опция – OutputFcn. В качестве параметра она принимает ссылку на функцию, которая и должна выполняться в процессе работы.

Эту функцию можно писать самому, однако есть несколько уже встроенных функций, к которым можно обращаться без объявления, а именно: 

• можно сразу строить графики всех компонент решения от времени, задавая функцию odeplot.
• строить график первых двух компонент решения посредством команды odephas2, что даже названием  намекает на фазовый портрет.
• с помощью команды odephas3 строится то же, только не для 2х а для 3х первых компонент решения.

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

Чтобы посмотреть, из чего состоят эти функции, воспользуйтесь командой edit, например введите editodeprint.

О том, по какому шаблону и как писать свои функции мы поговорим ниже. 

Давайте посмотрим, как задавать outputfcn в скрипте на конкретном примере.

В переменную опций через odeset запишем в параметр Outputfcn функцию odeplot через оператор собаку как при рассмотрении предыдущей опции. В качестве решаемой системы выберем уже известное нам по предыдущему ролику уравнение эйлера, даже название функции, заметьте, такое же.

Запустим скрипт.

Мы не использовали команду plot. Сама функция odeplot построила график. 

А теперь посмотрим работу odephas. Сменим odeplot на odephas2 и запустим скрипт.

Работает! Не сказал бы, что рисует красиво, однако вы в любом случае можете написать свои функции для отображения, которое Вас эстетически устраивает и уже скоро Вы узнаете, как их написать. 

Обратите внимание, что odeplot рисует СРАЗУ ВСЕ компоненты решения. А odephas2 рисует график из первых двух из ВСЕХ компонент решения. А если вам нужно только одно решение? Или из ста компонент решения 100-мерной системы вам нужно отображать и отслеживать только пару штук? 

На помощь приходит функция OutputSel. Она задаётся в паре с outputfcn. С помощью неё вы выбираете, какие компоненты решения попадут в выходную функцию. Вы задаёте вектор-строку номеров желаемых компонент. Более того, Вы можете формировать их порядок следования. 

В случае ситуации, показанной на рисунке, функция odeplot нарисует только 4-е и 2-е решения, а функция odephas2 по оси абсцисс будет откладывать x(4), а по оси ординат x(2), а не X(1) и X(2) соответственно, как это было бы без использования опции OutputSel. 

Посмотрим это в скрипте:

Добавим опцию outputSel и нарисуем odeplot.

Отобразился один график, знакомый нам по очертаниям по предыдущему ролику. 

Наконец переходим к “вишенке”.

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

В качестве аргументов этой функции следует задавать t – текущий аргумент. К нему можно обращаться из функции. Далее следует задать X – это вектор компонент решения. При этом, если опция refine больше единицы, то t будет представлять собой вектор-строку, чья длина будет равна значению refine, а X будет матрицей, где количество строк будет равно количеству компонент решения, а количество столбцов – значению опции refine, как показано на рисунке. 

Третьим аргументом является flag. Перед началом расчета, в этот аргумент передаётся значение ‘init’, после каждой итерации вычисления передаётся значение [] пустое, а после завершения расчета передаётся флаг ‘done’. Таким образом имеется возможность настроить функцию так, чтобы что-то происходило при инициализации, что-то после каждой итерации, а что-то, соответственно, после всего вычисления. Механику с флагами можно реализовать, например, через switch-case, но это уже кому как нравится. 

Функция должна возвращать переменную статуса. Её значение должно быть либо 0, тогда вычисления продолжатся, либо 1, тогда работа решателя останавливается.  

Например, функция может иметь такой вид:

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

Давайте рассмотрим простой вариант реализации кастомной outputfcn в живом скрипте.

Вот так выглядит функция. Перед началом вычисления она должна написать «Начало», в процессе должна писать выражение, заданное здесь, а в конце должна написать «Конец». 

Запустим.

Как хотели – так и получили. 

Функции ode и их опции можно, очевидно, использовать и в app designer. Рассмотрим такой пример.

В процессе расчета достраивается график. Расчет можно остановить, изменить начальные условия и запустить ещё несколько. Другой кнопкой можно очистить поле графика. Этот несложный пример выложен в открытом доступе на нашем репозитории.

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

Видеообзор по теме решения систем Д/У описанный в публикации доступен по ссылке.

Предыдущие публикации по теме: 

Часть 1. Решение систем обыкновенных дифференциальных уравнений в среде MATLAB

Часть 2. Решение систем обыкновенных дифференциальных уравнений в среде MATLAB. Опции решателей odeXY