Решение жестких систем дифференциальных уравнений

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

при на отрезке [0; 500].

Создадим вектор-функцию правых частей и вектор начальных условий:

Применим функциюRkadapt: . Система МС откажется выдать результат. Подсказка сообщает, что слишком много шагов интегрирования. Попробуем уменьшать отрезок интегрирования. Отрезки длиной 400, 300, 100, 50 дают такой же результат. Получим результат на отрезке [0; 20]:. Для вычислений компьютеру потребовалось значительное время. Построим графики компонент решения (рис. 41). Видим, что первая и третья компоненты решения резко возрастают до больших значений. Проверим систему на жесткость. Для этого создадим матрицу Якоби функцииF(x,y), т.е. матрицу из частных производных по компонентам y, пред дифференцированием скобки рекомендуется раскрыть, результат должен быть функцией от x и от y:

.

Найдем собственные числа этой матрицы в начальной точке решения . Так как все собственные числа отрицательны, то первое условие жесткости выполнено. Найдем число жесткости 0.161:0.001294 = 124.42. Это число не достаточно велико, чтобы признать систему жесткой. Посмотрим число жесткости в точке решения, где компоненты решения принимают большие значения:

1504:0.331 = 4559. Видим, что число жесткости очень большое. Следовательно, систему дифференциальных уравнений следует решать как жесткую систему.

Для решения жестких систем дифференциальных уравнений в МС14 предусмотрены несколько функций: AdamsBDF, BDF, Radau, Stiffb, Stiffr. Первая из них автоматически определяет, является ли система жесткой или нет, и выбирает соответствующий метод численного интегрирования (метод Адамса или формулу обратного дифференцирования). Остальные могут решать и нежесткие системы, но в этом случае объем вычислительной работы будет необоснованно завышен. Функции Stiffb, Stiffr поддерживаются значительно более ранними версиями МС, в то время как остальные появились в более позних версиях этой системы.

Рассмотрим использование функции AdamsBDF. Первые пять аргументов у нее такие же, как и у ранее используемых функций: вектор начальных условий, левая граница отрезка интегрирования, правая граница, число выводимых узлов интегрирования, имя функции в правой части системы уравнений. Шестой и седьмой аргументы являются необязательными. Шестой аргумент – имя матрицы Якоби (производные по компонентам неизвестной вектор-функции), эта матрица создается как функция двух аргументов. Седьмой аргумент – требуемая точность решения. Если седьмой аргумент больше или равен , то его значение на решение не влияет.

Найдем решение исходной системы дифференциальных уравнений . Компьютер практически моментально выполняет эту операцию. Посмотрим, с какой точностью МС выполнил интегрирование системы дифференциальных уравнений. Для этого добавим в функциюAdamsBDF необязательные аргументы. Матрица Якоби уже создана ранее и имеет имя A1. Отметим, что первым ее аргументом должно служить имя независимой переменной, как и было сделано ранее. Там где эта матрица создавалась, этот аргумент не требовался, матрица от него не зависит. Но она создавалась с учетом дальнейшего использования, когда этот аргумент необходим. Точность установим . Итак

Если проанализировать расхождение полученных результатов, то наибольшее различие наблюдается в узле с номером 608, т.е. при для первой компоненты решения. Это отличие составляет 26.654. В абсолютных единицах такое различие кажется большим, но если учесть что первая компонента решения в этой точке равна 83710, то ошибка составляет 0.03%, что вполне допустимо.

Построим графики решений (рис. 42). Число выводимых узлов интегрирования взято большим, тысяча, так как в противном случае на графиках могут быть пропущены пики возрастания функций. На рисунке графики второй и третьей компоненты практически не видны. Желающие могут построить график каждой компоненты решения отдельно, а также фазовые портреты решения на плоскостяхи.

Функция BDF имеет то же набор аргументов, что и функция AdamsBDF, Функция Radau имеет тот же набор обязательных аргументов и три необязательных, разбирать которые мы здесь не будем. Точность во всех трех этих функциях регулируется последним необязательным аргументом.

Функции Stiffb, Stiffr имеют шесть аргументов, причем все являются обязательными. Первые пять аргументов такие же, как и в описанных ранее функциях. Шестой аргумент – это имя матрицы Якоби функции в правой части системы дифференциальных уравнений, но частные производные берутся по всем аргументам, включая и независимую переменную. В этой матрице столбцов на один больше, чем строк.

Для системы уравнений, приведенной выше, она запишется так

.

Решение получим оператором или. Можно сравнить матрицыZ3 и Z4 с матрицей Z2. Для функций Stiffb, Stiffr точность результата регулируется системной переменной TOL.

Функция Stiffb использует метод Булирша–Штёра, Stiffr – метод Розенброка.