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

В данной публикации речь пойдет об опциях. Мы рассмотрим их синтаксис, рассмотрим опции уточнения, опции точности, научимся контролировать шаг и, что самое главное для вычислительно емких задач, поразмышляем относительно их обдуманного использования.

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

Как задавать опции?

Освежим в памяти синтаксис функций ode. Эти решатели можно настраивать. Для этого нужно в качестве последнего, четвертого аргумента функции ode подать переменную, содержащую опции:

Переменная эта создаётся с помощью другой функции odeset. В качестве аргументов ей следует подавать через запятую пары: название одной опции и её значение, название другой опции и её значение и так далее.

Теперь посмотрим на функцию odeset и создаваемую переменную в живом скрипте:

Как видно, переменная опций есть ни что иное как структура MATLAB. За что отвечают опции Reltol и Stats мы поговорим ниже. Всего имеется 22 опции. Некоторые из них требуются или могут быть заданы только для определённых решателей. Об опциях "общего назначения", подходящих ко всем решателям, мы поговорим в этом посте. 

Итак, с помощью функции odeset также возможно объединение двух переменных опций. Обе переменные можно объединить в одну переменную опций с заменой. Функции odeset следует подать 2 переменные опций. Результатом будет являться такая переменная опций, состоящая из опций второй переменной, а также из тех опций первой переменной, которые не были заданы во второй переменной опций. Так и осуществляется объединение с заменой.

Посмотрим пример:

Первая опция задаёт RelTol и Refine. 

Вторая опция задаёт Abstol и Refine с другим значением. 

Результатом объединения этих опций является новая переменная опций, которая, очевидно, содержит RelTol из первой переменной и Abstol из второй переменной, а вот опция Refine имеет значение второй переменной, ибо вторая переменная в данном случае является заменяющей.

Наконец-то приступаем непосредственно к рассмотрению опций. Поговорим о точностных опциях. Как и какую точность задать?

Решатели MATLAB различают 2 типа ошибок интегрирования:

1. относительную
2. абсолютную

Абсолютная ошибка AbsTol - модуль ошибки интегрирования.

Относительная ошибка RelTol - модуль отношения абсолютной ошибки к текущему значению.

По умолчанию, значение RelTol = 10^-3 и AbsTol=10^-6.

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

В процессе вычисления на всём интервале для каждой i-й компоненты решения справедливо выражение.

Шаг решателя подстраивается именно под эту величину. Однако можно заметить, что ошибка вычисляется для каждой компоненты отдельно, и для некоторых задач имеет смысл отслеживать норму ошибки. Чтобы это сделать, следует включить опцию NormControl (odeset ('NormControl' , 'on’)). По умолчанию, она выключена. Если её включить, то ошибка отслеживается многокомпонентно через норму, и выражение, по которому отслеживается ошибка вычисления, имеет вид, показанный на рисунке.

Здесь n –размерность решения.

С точностными опциями разобрались, двигаемся дальше.

Известно, что численные методы позволяют получить решения в узлах расчетной сетки. Однако может так случиться, что шаг интегрирования слишком большой, и, например, при визуализации хотелось бы иметь ещё больше точек на графике. В MATLAB существует возможность определить значение решения в промежуточных точках. Для этого используются так называемые "Формулы непрерывного расширения" (continuous extensionformulas), которые позволяют гладко интерполировать значения решений дифференциальных уравнений на промежуточные точки, причем точность значений в этих точках не ниже точности численного метода, заданного RelTol и AbsTol. 

Для подобного уточнения существует опция 'Refine’. В качестве параметра задаётся неотрицательное число (на сколько промежутков будет разделён интервал между узлами сетки).

Уже ранее показанная опция 'stats', как нетрудно догадаться из названия, показывает статистику решения. 

Перейдем в живой скрипт. Возьмём в качестве примера уравнение Эйлера.

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

Зададим начальные условия, а также зададим 2 переменных опций: в первой Refine будет 1, во второй – 100

Получим оба решения и заметим, что количество вызовов функций одинаково, что неудивительно.

Визуализируем координату x от времени для обоих решений. Сглаживание налицо.

Сейчас мы свернем немного в сторону от рассмотрения опций.

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

Но прежде следует осветить одну деталь. В предыдущие разы мы присваивали переменным t и X результат выполнения функции решателя ode. Если присвоить одной переменной значение функции решателя, то эта переменная будет структурой, содержащей в себе информацию о решении.

Для определения значений решения в желаемых точках на помощь приходит  функция deval со следующим синтаксисом:

Часть 1

Часть 2

Часть 3

Заметьте, что функция называется deval а не, скажем, odeval, потому, что она также работает с результатами другого типа решателей - dde (delay differentialequation) или уравнения с запаздыванием. 

Теперь же вернемся к рассмотрению опций, а именно опций контроля шага.

Большинство решателей MATLAB определяют размер шага автоматически, исходя из заданных точностных опций. Тем не менее, имеется возможность "насильно" ограничить шаг решателя, задавая максимальный шаг с помощью опции 'MaxStep’.

По умолчанию MaxStep есть одна десятая от всего временного интервала: 0.1*abs(t0-tend).

Также можно задать начальный шаг. Если его специально не указывать, то решатель выбирает начальный шаг, исходя из наклона компонент решения в начальный момент времени. Однако, если, например, наклон всех компонент решения нулевой, то решатель выберет слишком большой размер шага. Для избегания подобной ситуации и следует задавать начальный шаг с помощью опции initialStep.

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

Видно, что уже на 3-м шаге интегрирования шаг решателя "упёрся" в заданное нами значение, что свидетельствует о том, что заданную точность можно было бы достичь и большим шагом. И вообще, нужно ли искусственно управлять шагом решателя если есть возможность задавать точность и уточнять решение? Рассмотрим это подробнее. 

Сравним все 3 подхода на решателе ode23. Возьмем за основу уже использовавшийся ранее дифур. 

• Для первого случая снизим максимальный размер шага.
• Для второго случая выполним уточнение через Refine.
• Для третьего случая увеличим точность через AbsTol и RelTol.

Параметры подберем таким образом, чтобы длины выходных векторов были близки друг к другу. Будем засекать время вычисления с помощью команд tic-toc.

Как можно заметить, контроль за максимальным шагом не даёт видимых преимуществ для большинства задач.

Если вам требуется уточнить решения для, например, постройки графиков, то используйте refine, а не уменьшайте шаг. Вам это не даст преимуществ, и лишь увеличится время расчета. 

Если же Вам требуется высокая точность, то используйте точностные функции abstol и reltol. Искусственное уменьшение шага не влияет на точность вычисления никоим образом, это частое заблуждение тех, кто пользуется решателями ode. Не заблуждайтесь! Ограничивайте шаг только в самом крайнем случае, когда решатель может пропустить какую-то особенность решения. Особенно это касается решения задач умеренной жесткости нежёсткими решателями, но эта тема выходит за пределы этой публикации.

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