Мульти-неявные методы со второй производной для жестких систем дифференциальных уравнений Текст научной статьи по специальности «Математика»

© Васильев Е.И., Васильева Т. А., Киселева М.Н., 2012

УДК 519.62 ББК 22.19

МУЛЬТИ-НЕЯВНЫЕ МЕТОДЫ СО ВТОРОЙ ПРОИЗВОДНОЙ ДЛЯ ЖЕСТКИХ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Е.И. Васильев, Т.А. Васильева, М.Н. Киселева

Представлено новое семейство А-устойчивых разностных схем для численного решения жестких систем обыкновенных дифференциальных уравнений. Основным отличием семейства является мульти-неявность с использованием вторых производных искомого решения. Подробно рассмотрены схемы 4-го, 6-го и 8-го порядков аппроксимации. Приведены результаты тестирования разностных схем на задаче Хилла и тесте Крейса. Представлена процедура автоматического выбора шага интегрирования с контролем точности, основанная на комбинации двух методов. Показано значительное преимущество в эффективности схем высокого порядка с автоматическим выбором шага интегрирования.

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

Введение

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

При решении жестких систем возникает ряд проблем [4-11], связанных с устойчивостью, точностью и вычислительной сложностью применяемых методов. Наиболее полные обзоры по численным методам для жестких систем содержатся в [4; 10]. Центральной в теории численных методов решения жестких систем дифференциальных уравнений является проблема устойчивости. Она заключается в сложности одновременного достижения высокого порядка аппроксимации и абсолютной устойчивости разностной схемы. Хорошо известно, что явные методы плохо применимы для жестких систем, так как не являются А-устойчивыми, а для классических неявных методов существует так называемый порядковый барьер Далквиста [6; 7; 9] - неявные линейные многошаговые методы выше второго порядка не могут быть А-устойчивыми. В неявных методах поиск решения на каждом шаге интегрирования осуществляется итерационно. Для жестких систем, чтобы не нарушить устойчивость, обычно применяется метод Ньютона, в котором требуется вычислять матрицу производных правой части системы по искомому решению. В неявных методах Рунге-Кутта [8] метод Ньютона на каждом шаге интегрирования требуется применять несколько раз для определения неизвестных вспомогательных значений правой части. По сути эти методы являются несколько раз неявными с многократным вычислением матриц производных при реализации метода Ньютона. В принципе эти матрицы производных можно бы использовать и на заключительном этапе метода с целью повышения порядка аппроксимации. Таким образом, имеет смысл рассмотреть в прямом

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

В данной работе реализуется подобный подход построения семейства мульти-неявных методов со вторыми производными решения. Вторые производные приводят к включению в разностную схему матрицы Якоби правой части. Проведен анализ аппроксимации и устойчивости некоторых методов семейства. Представлен алгоритм автоматического выбора шага интегрирования с контролем точности. Предложенные методы были протестированы на задаче Хилла [1] и тесте Крейса [11]. Сравнение результатов показало значительное преимущество в эффективности схем высокого порядка с автоматическим выбором шага интегрирования.

1. Мульти-неявные методы со второй производной

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

и(7) = f(и), и(0)= и0, t> 0, (1)

где и = {и1(7), и2(7), ... , ир(7)}т и f(и) являются векторными функциями размерности р. Введем сетку wt = {7 = тп, п = 0, 1, 2...} с постоянным шагом т> 0. Функции, определенные на сетке, и приближенные решения обозначим как уп = у(7п), / = f (уп).

Для численного решения задачи (1) предлагается т раз неявный разностный метод

Уп+ *~ Уп+ к-1 = ет Ш + 7Ьуп+,)/п+,,к = 1,2,...,т, (2)

7 ,= 0

где у1 и/р - векторы. Е и J- единичная матрица и матрица Якоби размерархр:

/ = /(Уг ), Jг = 7/ (У, X 0' = n, п +1,..., п + т).

^ и

При известном уп = у(7п) разностная схема (2) представляет собой систему взаимосвязанных нелинейных разностных уравнений относительно неизвестных уп+1, Уп+2, ... Уп+т в т последующих точках 7п+1, 7п+2, ... 7п+т . Таким образом, мы имеем т раз неявную схему.

Для нахождения коэффициентов разностной схемы (2) запишем ошибку аппроксимации, используя свойство решения системы (1) и у= /ии у= /и/:

и - и т

ук = - п+к п+к-1 + е (аки%, + 7Ъ'ки%,), к = 1,2,...,т . (3)

7 ,= 0

Разложение функций, первых и вторых производных в (3) в ряд по степеням тс центром в точке 7п позволяет выписать систему линейных уравнений на коэффициенты, которые обеспечивают

уп = 0(72т+2), к = 1,2,...,т .

Ниже приведены коэффициенты для случаев т = 1, 2 и 3 в матричной форме записи.

При т = 1:

ж1 1ц ж1 11

(а)=Ё 1(Ъ‘)-^ - 721

При m = 2:

При m = 3:

К )=ттт:

і ж01

240 иіі

128

128

11ц

і0і!

1 ж13

240 и3

40

40

- 3ц

- 13І

(4 ) =

18144

ж6893 8451 2403 397ц

1243 8829 8829 2431,

и397 2403 8451 б893ш

(bk ) =

30240

Й283 - 7б59 - 2421 - 163ц

193 3051 - 3051 9 -

Шб3 2421 7б59 - 1283 ш

1

1

Разностные схемы с этими коэффициентами имеют соответственно 4-й, б-й и 8-й порядок аппроксимации. Две последние впервые были исследованы в работах [3; 12]. Далее будем их именовать MISD4, MISD6 и MISD8 соответственно (Multi-Implicit Second Derivative).

Исследуем схему (2) на устойчивость. Разностная схема называется абсолютно устойчивой, если численное решение модельного уравнения (1) с/(u) = lu (1 комплексное, Re1<0) асимптотически стремится к нулю с ростом t при любой величине шага интегрирования t Применяя схему (2) к модельному уравнению, находим отношение уг&к к yn , которое будет зависеть от произведения z = 1t

— = Rk (2)

у к V ✓ .

У п

Область в комплексной плоскости, в которой функция роста ^.г)|<1, называется областью устойчивости. Для А-устойчивости необходимо, чтобы область устойчивости включала в себя всю левую полуплоскость.

При т = 1

R1( z ) =

При m = 2 функция R(z) найдена в [12]

R2( z) =

1+ z +

104 z 2 + 240

ті z 3 +

1- z+

104 z - -1 z + -1 z

240 10 ^ 90

(4)

При m = 3 функция R(z) найдена в [3]

R3( z ) =

rz 4 +

z+

rz 4-

z 5 + z 6

560 560

1

Во всех приведенных случаях числитель и знаменатель отличаются только знаком при нечетных степенях. Это приводит к тому, что при чисто мнимых значениях 2 числитель и знаменатель являются комплексно сопряженными, то есть модули равны. Следовательно, границей области устойчивости является мнимая ось комплексной плоскости. Также видно, что для отрицательных действительных 2 во всех случаях ^.г)| < 1. Таким образом, область устойчивости в рассмотренных случаях совпадает с левой полуплоскостью комплексной плоскости 2. Следовательно, разностная схема (2) при т J 3 А-устойчива.

Заметим, что это свойство выполняется и для произвольного т, так как схема (2) и по структуре и по найденным коэффициентам обладает некоторой симметрией. Она инвариантна при замене 7 ® -7, т® -т, что приводит к взаимной симметрии областей устойчивости и неустойчивости.

2. Тестирование мульти-неявных разностных схем

На рисунке 1 представлены результаты тестирования рассматриваемых методов на задаче Хилла [1], которое было выполнено в работе [3]. Выполнялось интегрирование системы из четырех дифференциальных уравнений первого порядка на небольшом участке по времени 7 О [0, 2] с различными шагами интегрирования т и вычислялась норма разности численных решений с конт-6ТёШ иI (оТЧ и! ) баз а еа а ёТ Г и е I ТI а о 7 = 2. Контрольное решение предварительно рассчитывалось с максимально возможной точностью.

0 -1 -2 -3 І8(г) -4 0 1 2 3 4

а б

Рис. 1. Зависимость погрешности є численного решения: а) от величины шага интегрирования т б) от затрат времени СРи. Z0 - некоторый масштаб затрат времени СРи

На рисунке 1а в логарифмических осях хорошо видно, что точность методов соответствует теоретической, за исключением начального и конечного участков для методов высокого порядка. Хаотическое поведение погрешности в районе е = 10-14 объясняется пределом машинной точности. Данные расчеты проводились со стандартной двойной точностью формата вещественных чисел с мантиссой 52 бит. Нарушение линейности на начальном участке при относительно больших шагах т объясняется вкладом погрешностей более высокого порядка.

Нахождение решения на каждом шаге интегрирования в MISD-схемах, как и во всех неявных схемах, осуществляется итерационно. Для жестких систем, чтобы не нарушить устойчивость, обычно применяется метод Ньютона, в котором требуется вычислять матрицу производных правой части по искомому решению. Но сами матрицы производных уже содержатся в правой части MISD-схем, поэтому использование метода Ньютона приводит лишь к малому дополнительному увеличению вычислительных затрат. Практика показала, что вычислять производные от матриц Jn+г. в правой части формул (2) нет необходимости из-за присутствия перед ними малого коэффициента. Даже в сильно жесткой задаче химической кинетики [2] это не влияло на сходимость метода Ньютона, для которой обычно достаточно максимум 3 итераций. В этом случае основные вычислительные затраты на каждом шаге интегрирования для т-неявного метода включают: 1) 3т раз вычислений матрицы производных размера р хр; 2) 3т раз перемножений матриц размера р х р ; 3) 3 раза решений системы линейных уравнений размера тр г тр. Затраты первого и второго вида пропорциональны т, затраты третьего вида пропорциональны т3. Линейные по т затраты преобладают в системах уравнений невысокой размерности с громоздкими нелинейными выражениями в правой части. Кубические по т затраты преобладают в системах уравнений высокой размерности с разреженной матрицей производных.

Хотя в прикладных задачах обычно имеет место первый случай, оценим эффективность m-неявных методов между собой для второго случая. В этом случае вычислительные затраты на каждом шаге интегрирования для MISD8, MISD6 и MISD4 будут относиться как 27 : 8 : 1. Поэтому для сравнительной оценки затрат в логарифмической шкале при условии одинаковой точности расчета нужно сместить графики на рисунке 1а вправо на величину lg(27) для MISD8 и на lg(8) для MISD6. Смещенные графики изображены на рисунке 16. Видно, что метод 6-го порядка в широком диапазоне точности эффективнее метода 4-го порядка. Метод 8-го порядка эффективнее метода 6го порядка лишь при необходимости очень точных вычислений при допустимой погрешности порядка 10-12 и менее. Для получения максимальной отдачи от метода MISD8 требуется аппаратно или программно обеспечить более точную машинную арифметику с вещественными числами.

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

3. Тест Крейса

На рисунках 2-3 представлены результаты тестирования метода MISD6 на жесткой системе задачи Крейса [11], которое было выполнено в работе [12]. Тест Крейса представляет собой неавтономную систему двух дифференциальных уравнений:

— ж e 1 пи

—u(t)= E(t)£ JE"ЧОг U(t), (5)

где u(t) = (u1(t), U2(t))T,

acos t - sin t

E (t )=L t t

item t cos t

Начальные условия u(0) = (-0,7; 0,7)T; e = 0,05. При численном решении система преобразовывалась к автономной системе трех дифференциальных уравнений.

Графики решений (см. рис. 2) имеют две типичные зоны поведения решений жесткой системы. Первая зона включает быстрое изменение функции ux(t), когда она переходит к равновесному значению. Затем эта функция медленно меняется вблизи равновесного значения. В этом примере равновесное значение зависит от времени t.

Приближенное решение системы (5) обозначим y(t) = (yx(t), y2(t)). Расчеты проводились с постоянным шагом. Отрезок интегрирования [0; 3] разбивался на N равных частей с t = 3/N. Оценивалось отклонение v(t) = 1 (t) — ux(t)| от контрольного (точного) решения на всем промежутке интегрирования. Рисунок 3 демонстрирует поведение нормированной на t6 ошибки для MISD6-метода и для метода Гира 6-го порядка (см.: [9]). Обратим внимание, что для метода Гира шкала по вертикали в 100 000 раз крупнее.

Из сравнения этих графиков можно заключить, что константа погрешности у метода MISD6 в 100 000 раз меньше, и для достижения такой же точности в методе Гира необходимо использовать шаг интегрирования в 100 0001/6 » 7 раз мельче.

С другой стороны, вычислительные затраты на один шаг по времени в методе MISD6 несколько выше. Основные затраты в этой задаче уходят на вычисление матриц производных: 6 матриц в методе MISD6 против 3 матриц в методе Гира. С учетом дополнительных накладных расходов на перемножение матриц и решение систем линейных уравнений метод MISD6 требует в среднем в 3 раза больше времени CPU на один шаг интегрирования по сравнению с методом Гира. Следовательно, при равномерном шаге интегрирования эффективность MISD6 в 2 раза выше, то есть для расчетов с одинаковой точностью требуется приблизительно в 2 раза меньше времени CPU. Однако с учетом сложности программирования можно признать эффективность обоих методов приблизительно равной.

0< e <<1; 0J t J 2p.

Щей

О 1 2 т 3

Рис. 2. Решение задачи Крейса

/N=300 ^N=6 00 Л^=1200

N=2400

О 1 2 I 3 О 1 2 f 3

а б

Рис. 3. Нормированные погрешности метода MISD6 (а) и метода Гира (б) на равномерной сетке

4. Выбор шага интегрирования с контролем точности

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

Погрешностью метода является разность между точным и приближенным решениями исходного уравнения (1). Но точное решение не всегда известно. Рассмотрим две разностные схе-

мы с 6-м и 8-м порядком аппроксимации. Различие в два порядка аппроксимации позволяет принять приближенное решение 8-го порядка в качестве точного решения уравнения (1) по сравнению с приближенным решением 6-го порядка. Оцениваем разность между решениями этих разностных схем и по ней контролируем шаг интегрирования. Для контроля точности из (2) используем разностное уравнение, получаемое суммированием всех разностных уравнений. Оно связывает значения в точках t и ^ .

п п+т

Пусть у = {уп, уп+1, Уп+2, Уп+3} - решение, полученное по схеме MISD8. Подставим его в правую часть суммарного уравнения схемы MISD6 вместе с уже найденными значениями Ау) и матрицами производных и вычислим vn+2:

V - У

^+"2^ = F(Уn , Уп+^ Уп+ 2 ) . (6)

Вычисляем разность В = || уп+2 - vn+21|. Из (6) с учетом 6-го порядка аппроксимации схемы MISD6 следует

^+22/п+2 = F(Уп,Уп+1,Уп+2)- Уп±^ = 0(г6).

В результате:

1^п+2- Уп+2 11= В =2С 7. (7)

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

Предположим, что интегрирование должно выполняться на отрезке [0; Т] с допустимой ошибкой е для всего интервала. Тогда допустимую ошибку на шаге 2 т следует принять за 2ет/Т. Поэтому новый (следующий) шаг интегрирования выбираем из условия

7 2t 2Ct = e-

Т

где используется С с текущего шага. Подставляя выражение С через В из (7), получим связь между величиной текущего и следующего (нового) шага интегрирования:

t \2te

= б-----. (8)

t Чтв ’

Таким образом, алгоритм автоматического выбора шага состоит из следующих стадий:

1) Vn, yn+1, yn+2, Vn+3 находятся по схеме 8-го порядка (MISD8);

2) из формулы (б) вычисляем vn+2 и затем B = || yn+2 - vn+21|;

3) следующий шаг интегрирования tnew вычисляется по формуле (8).

Чтобы избежать слишком резкого увеличения или уменьшения шага интегрирования, разумно накладывать ограничения на рост или уменьшение шага. Обычно используются:

1 < t new < 2

2 ~Г .

Вышеизложенный алгоритм изменения шага интегрирования без труда обобщается и на

другие пары MISD-методов. Например, для расчета по схемам 8-го или б-го порядков с

контролем шага по методу 4-го порядка. Заметим, что затраты времени на контроль точно-

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

функций или матриц производных.

а б

Рис. 5. Погрешность MISD6 с контролем точности: а) е = 2х10-9; б) е = 3х10-11

На рисунке 5 показаны результаты расчетов с постоянным и переменным шагом интегрирования для задачи Крейса. Расчеты проводились по схеме 6-го порядка, контроль и изменение шага осуществлялись по схеме 4-го порядка. Изображены абсолютные погрешности первой компоненты решения. Непрерывная линия есть результаты с постоянным шагом т= 3/Ж Линии с точками есть результаты с переменным шагом, около них указано количество шагов, которое потребовалось на весь отрезок интегрирования. Величина шага рассчитывалась по формуле

Погрешность £ задавалась как максимальная погрешность, возникающая при постоянном шаге интегрирования, она указана в подписи к рисунку 5. Видно, что в среднем процедура выбора шага уменьшает размер расчетной сетки (и затраты машинного времени) более, чем в 4 раза при сохранении требуемой точности.

Заметим, что метод Гира, как и большинство линейных многошаговых методов высокого порядка, имеет весьма ограниченный ресурс изменения величины шага интегрирования. Контроль точности для него организовать несложно, но при изменении шага всего лишь на несколько процентов метод Гира теряет свои свойства устойчивости. MISD-методы таким недостатком не обладают, так как не используют информацию о предыстории решения и, в принципе, могут изменять шаг произвольным образом.

На рисунке 6 проводится сравнение процедур выбора шага с контролем точности для трех различных пар MISD-методов. В первой паре расчет проводился по схеме 6-го порядка с контролем шага по схеме 4-го порядка, во второй паре расчет проводился по схеме 8-го порядка с контролем по схеме 4-го порядка и в третьей паре расчет проводился по схеме 8-го порядка с контролем по схеме 6-го порядка. Все расчеты проводились с постоянным шагом т = 0,08, новый шаг %пем! вычислялся, но не применялся.

О 0.4 0.8 1.2 1.6 і 2.0

Рис. 6. Коэффициент изменения шага при различных комбинациях методов, т = 0,08; є = 3 х 10-6

Таким образом, мы можем сравнить максимально возможное изменение шага интегрирования в течение расчета для различных пар. Видно, что коэффициент допустимого изменения шага зависит как от порядка точности расчетной, так и от порядка контролирующей схемы. Однако порядок точности контролирующей схемы влияет гораздо больше. На большей части отрезка интегрирования контроль по схеме 6-го порядка позволяет увеличить шаг в 8 раз, в то время как контроль по схеме 4-го порядка - не более чем в 1,5 раза. Учитывая ранее сделанные оценки затрат MISD-схем на одном шаге интегрирования, видим, что пара (8 - 6) по вычислительным затратам приблизительно в 3 раза эффективнее пары (6 - 4).

Заключение

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

СПИСОК ЛИТЕРА ТУРЫ

1. Батхин, А. Б. Порождающие плоские периодические орбиты задачи Хилла / А. Б. Батхин // Препринты Ин-та приклад. мат. им. М. В. Келдыша РАН. - М., 2010. - N° 47. - 24 с. - Режим доступа: http://library.keldysh.ru /preprmt.asp?id=2010-47.

2. Васильев, Е. И. Влияние инертных частиц на структуру стационарной волны детонации в гремучей смеси (2Н2 + 02) / Е. И. Васильев, А. А. Васильченко, Т. А. Васильева // Вестн. ВолГУ Сер. 1, Математика. Физика. - 2002. - Вып. 7. - С. 54-65.

3. Киселева, М. Н. Неявная А-устойчивая схема восьмого порядка для жестких систем дифференциальных уравнений / М. Н. Киселева // Выпускная квалификационная работа по направлению подготовки «Прикладная математика и информатика». - Волгоград, 2011. - 34 с.

4. Ракитский, Ю. В. Численные методы для решения жестких систем / Ю. В. Ракитский, С. М. Устинов, И. Г. Черноруцкий. - М. : Наука, 1979. - 208 с.

5. Самарский, А. А. Численные методы / А. А. Самарский, А. В. Гулин. - М. : Наука, 1989. - 432 с.

6. Dahlquist, G. Error analysis of a class of methods for stiff non-linear initial value problems / G. Dahlquist // Numerical Analysis. Lecture Notes in Mathematics 506. - Berlin : Springer, 1975. - P 60-74.

7. Dahlquist, G. A Special Stability Problem for Linear Multistep Methods / G. Dahlquist // BIT. - 1963. -№3.- P. 27-43.

8. Dekker, K. Stability of Runge-Kutta methods for stiff nonlinear differential equations / K. Dekker, J.G. Verwer. - Amsterdam : Elsevier Science Publishers, 1984. - 318 p.

9. Gear, C. W. ODE methods for the solution of differential/algebraic systems / C. W. Gear, L. R. Petzold // SIAM J. Numer. Anal. - 1984. - № 21. - P. 716-728.

10. Hairer, E. Solving ordinary differential equations. II. Stiff and differential-algebraic problems / E. Hairer, G.Wanner. - Berlin : Springer-Verlag, 1991. - 614 p.

11. Kreiss, H. O. Difference methods for stiff ordinary differential equations / H. O. Kreiss // SIAM J. Numer. Anal. - 1978. - № 21. - P. 21-58.

12. Vasilev, E. High order implicit method for ODEs stiff systems / E. Vasilev, T. Vasilyeva // Korean Journal of Computational & Applied Mathematics. - 2001. - V 8, № 1. - 165-180.

MULTI-IMPLICIT SD-METHODS FOR STIFF SYSTEMS OF DIFFERENTIAL EQUATIONS

E.I. Vasilev, T.A. Vasilyeva, M.N. Kiseleva

The new set of absolutely stable difference schemes for a numerical solution of ODEs stiff systems is submitted. The main feature of the set is the multi-implicit finite differences with the second derivatives of the desired solution. The schemes with 4, 6 and 8 orders of accuracy in detail are considered. The testing results of difference schemes on the Hill’s problem and the Kreiss stiff system are shown. The procedure of automatic selection of an integration step with the control of accuracy is submitted. The considerable advantage in efficiency of the high order schemes with the automatic step selection is shown.

Key words: finite difference schemes, stiff systems, implicit methods, multi-implicit methods, methods with second derivative, accuracy control, step selection.