Алгоритм интегрирования с переменным числом стадий для решения умеренно жестких задач

Введение

При численном исследовании жестких задач вида

, , ,                                                   (1)

где y и f – гладкие вещественные N-мерные вектор функции; t – независимая переменная; все большее внимание привлекают явные методы [1–2]. Это связано с тем, что при применении L-устойчивых методов возникает проблема с обращением матрицы Якоби. В случае большой размерности системы дифференциальных уравнений время декомпозиции данной матрицы фактически определяет общие вычислительные затраты. В то же время явные методы не нуждаются в вычислении матрицы Якоби и, если жесткость задачи не слишком велика, то они будут предпочтительнее. В настоящее время применяются алгоритмы переменных порядка и шага, что приводит к существенному повышению эффективности расчетов. На участках установления нет смысла использовать численные формулы высокого порядка точности. Быстродействие можно повысить за счет применения методов низкого порядка, но с большими областями устойчивости. Поэтому дальнейшее повышение эффективности достигается за счет построения алгоритмов интегрирования не только с переменными порядком и шагом, но и с переменным числом стадий.

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

Явные методы типа Рунге-Кутта

Для численного решения жестких задач в [2] предлагается применять явные методы вида

,   ,                                     (2)

где k_i, 1≤i≤m – стадии метода; p_mi, α_ij и β_ij – коэффициенты, определяющие свойства точности и устойчивости схемы (2). Для упрощения выкладок далее рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений

,   ,   ,                                                  (3)

для решения которой применим методы вида

,   ,   ,                        (4)

где k_i=hf(y_n,i-1), 1≤i≤m, y_n,0=y_n. Все полученные ниже результаты можно использовать для неавтономных задач, если в (2) положить

,  ,  .                                                       (5)

Условия порядка

Ниже потребуется матрица B_m с элементами b_ij [1–2]

,   ,

,  ,  ,                                                                  (6)

,  ,  ,

где β_ij – коэффициенты схемы (2) или (4). Устойчивость одношаговых методов обычно исследуется на линейном скалярном уравнении

,  ,                                                                 (7)

с комплексным λ, Re(λ)<0. Применяя вторую формулу (4) к (7), получим

, , , , .                (8)

Обозначая C_m=(c_m1, …, c_mm)^T и P_m=(p_m1, …, p_mm)^T, последнее соотношение (8) можно переписать в виде

,                                                                                       (9)

где элементы матрицы B_m определены соотношениями (6). Для промежуточных численных схем (4) имеем

,  , ,  ,      .        (10)

Используя обозначения β_k=(β_k+1,1, …, β_k+1,k)^T и c_k=(c_k1, …, c_kk)^T, получим, что коэффициенты β_ij промежуточных схем (4) и коэффициенты соответствующих многочленов устойчивости связаны соотношениями

                                                                (11)

Заметим, что из сравнения (6) и (10) следует, что b_ki=c_i-1,k-1, то есть элементы (k+1)-го столбца матрицы B_m совпадают с коэффициентами многочлена устойчивости Q_k(z). Отсюда следует, что если заданы коэффициенты многочленов устойчивости основной и промежуточной численных схем, то коэффициенты методов (4) однозначно определяются из линейных систем (9) и (11) с верхними треугольными матрицами B_i, 1≤i≤m.

Разлагая точное и приближенное решения в ряды Тейлора по степеням h, можно записать

,

                                        (12)

,

где элементарные дифференциалы вычислены на точном y(t_n) и приближенном y_n решениях, соответственно. Из сравнения соотношений (12) при условии y_n=y(t_n) видно, что численная формула (4) будет иметь первый порядок точности, если

.

Отсюда следует, что для построения m-стадийных методов первого порядка точности в линейной системе (9) следует положить c_m1=1.

Для того чтобы воспользоваться соотношениями (9) и (11) требуются коэффициенты многочленов устойчивости. Функция устойчивости явного m-стадийного метода типа Рунге-Кутта представляет собой полином степени m. Через заданные коэффициенты схемы (4) можно вычислить коэффициенты этого многочлена по формуле (9).

Многочлены устойчивости

Пусть заданы два целых числа k и m, k≤m. Рассмотрим многочлены вида

,                                                          (13)

где коэффициенты c_i, 1≤i≤k заданы, а c_i, k+1≤i≤m − свободные. Обычно c_i, 1≤i≤k определяются из требований аппроксимации. Поэтому для определенности ниже будем предполагать, что c_i=1/i!, 1≤i≤k.

Обозначим экстремальные точки (13) через x₁, … ,x_m−1, причем x₁>x₂>…>x_m−1. Неизвестные коэффициенты c_i, k+1≤i≤m определим из условия, чтобы многочлен (13) в экстремальных точках x_i, k≤i≤m−1 принимал заданные значения, то есть

, ,                                                               (14)

где F(x) есть некоторая заданная функция, F_i=F(x_i). Для этого на x_i, k≤i≤m−1 и c_j, k+1≤j≤m рассмотрим следующую алгебраическую систему уравнений

, , , .                         (15)

Перепишем (15) в виде, удобном для расчетов на ЭВМ. Для этого обозначим через y, z, g и r векторы с координатами

, , ,

, ,                                                                 (16)

через E₁,…,E₅ − диагональные матрицы с элементами на диагонали вида

, , ,

,                          (17)

а через A – матрицу с элементами , . Заметим, что компоненты векторов (16), матриц (17) и  зависят от чисел  и , причем

,  ,  ,  ,  ,  .

Ниже аргументы опущены для упрощения записи. С использованием введенных обозначений задачу (15) можно записать в виде

,   .                                                                       (18)

Система (18) плохо обусловлена, что приводит к определенным трудностям при использовании для ее решения метода простой итерации. Для сходимости метода Ньютона требуются хорошие начальные условия, что в данном случае есть трудновыполнимая проблема. Если в (15) положить F_i=(−1)ⁱ, k≤i≤m−1, то есть поставить задачу нахождения полинома с максимальным размером интервала устойчивости, то вопрос о вычислении начального условия y⁰ решается достаточно просто с использованием значений экстремальных точек многочлена Чебышева, рассматриваемого на отрезке [−2m²,0], где m – есть степень (13). Их можно вычислить по формуле

,   .                                             (19)

Подставляя (19) в первую формулу (18), получим коэффициенты полинома Чебышева, для которого |Q_m1(x)|≤1 при xÎ[−2m²,0]. При любом  в качестве начальных условий можно взять (19) и, как показывают расчеты, имеется хорошая сходимость. Если же F_i≠(−1)ⁱ, k≤i≤m−1, то выбор начальных условий является, вообще говоря, искусством.

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

.                                                (20)

Ясно, что после определения стационарной точки (20) коэффициенты полинома устойчивости можно вычислить из первой системы (18). Заметим, что при использовании матрицы E₅ все собственные числа матрицы Якоби задачи (20) имеют отрицательные вещественные части, то есть задача (20) устойчивая. Из результатов расчетов следует, что (20) является жесткой задачей. Методы решения таких задач предполагают вычисление матрицы Якоби, что в случае (20) связано с трудностями. Поэтому для ее решения используем метод второго порядка точности с численным вычислением и замораживанием матрицы Якоби [3−4].

Метод первого порядка

В [5] описан алгоритм построения многочленов с заданными свойствами на промежутке [−1, 1], который позволяет вычислить коэффициенты многочленов устойчивости, соответствующих методам первого порядка, до степени m=27. Если параметры схемы (2) и коэффициенты полинома (13) связаны соотношениями (9) при условии c_m1=1, то решение (9) дает параметры m-стадийной схемы (2) первого порядка точности. При заданных β_ij, 2≤i≤m, 1≤j≤i−1 и c_ij, 2≤i≤m задача (13) имеет единственное решение в силу невырожденности матрицы B_m. При выборе β_ij можно руководствоваться различными соображениями. Если, например, набор формул (2) используется совместно с методом более высокого порядка точности, то и в схемах первого порядка естественно применять те же самые k_i.

Получим неравенства для контроля точности и устойчивости. Схемы первого порядка предполагается использовать на участке установления, где шаг ограничен устойчивостью, а не точностью. Контроль точности в этом случае носит вспомогательный характер, и для получения соответствующего неравенства будем использовать оценку локальной ошибки. В [2] показано, что из сравнения точного решения y(t_n+1) задачи (1) и приближенного решения y_n+1, полученного по схеме (2), до членов h² включительно при условии y_n=y(t_n) локальная ошибка схемы (2) первого порядка точности имеет вид

.                                             (21)

Учитывая (8), ее можно переписать так:

,                                                      (22)

где c_2m – заданный коэффициент при x² в (13).

Величину δ_n,1с помощью уже вычисленных стадий k_i, 1≤i≤m можно оценить многими способами. Будем поступать следующим образом. Введем обозначения:

, ,                                      (23)

где

, .                                                  (24)

Тогда, учитывая что

,  ,                      (25)

для контроля точности вычислений и при выборе шага интегрирования можно применять соответственно неравенства

,.                                                                      (26)

Вектор k₁ зависит от размера шага линейно. Поэтому с помощью первого неравенства (26) повторное вычисление решения сопровождается всего лишь одним дополнительным вычислением правой части дифференциальной задачи. В случае успешного шага второе неравенство (26) не приводит к увеличению вычислительных затрат, потому что вектор f(y_n+1) не применяется в следующем шаге интегрирования. Если второе неравенство (26) используется для контроля точности вычислений, то возврат при невыполнении точности будет дорогим, причем с ростом m затраты увеличиваются. Однако предварительный контроль A'_n позволяет в основном избежать повторных вычислений решения.

В случае применения численных формул (2) с переменным числом стадий при 2≤m≤M, где M – заданное целое число, можно поступать следующим образом. Введем обозначения:

,  ,                                                    (27)

где g''_m1вычисляются по первой формуле (24). Оценку ошибки определим по формулам

 , .                                    (28)

Тогда для контроля точности и при выборе шага интегрирования сразу для всего набора формул снова можно использоваться неравенства (26).

Теперь перейдем к получению неравенства для контроля устойчивости. Ниже будем предполагать выполненными соотношения

,  ,                                                                     (29)

где β₂₁, β₃₁, β₃₂ – параметры схемы (2).

Неравенство для контроля устойчивости получим следующим образом. Запишем разность (k₃ – k₂) в форме Тейлора с остаточным членом в форме Лагранжа первого рода:

.                                                        (30)

Учитывая вид разности (k₂ – k₁), оценку максимального собственного числа матрицы Якоби можно определить степенным методом, т.е. для контроля устойчивости методов (2) первого порядка будем применять неравенство

,                                                                                (31)

где

,                                                                (32)

а положительные постоянные γ_m,1, 3≤m≤M связаны с размером области устойчивости численных схем. Неравенство (31) можно использовать на каждом шаге для выбора эффективной численной схемы.

Таким образом, задача о построении явных методов Рунге-Кутта первого порядка точности с заданной областью устойчивости сводится к решению системы линейных алгебраических уравнений (9) с невырожденной матрицей B_m, где компоненты вектора c_m определяют размер и форму области устойчивости. Ниже с помощью построенных численных схем сформулируем алгоритм интегрирования, в котором допускаются расчеты, как по фиксированной схеме, так и с переменным числом стадий. Параметры методов можно хранить в виде двухмерного массива E с элементами e_ij вида

,  ,  ,  ,  ,  ,                     (33)

где M – максимальное число стадий в наборе методов. Величины γ_m,1, 3≤m≤M будем хранить в одномерном массиве.

Алгоритм с переменным числом стадий

Пусть имеется набор m-стадийных методов, длины интервалов устойчивости которых равны γ_m,1, 3≤m≤M, M ≥ 3. Ниже будем полагать, что при первом входе m=3. Алгоритм интегрирования с переменным числом стадий на основе методов первого порядка записывается следующим образом.

Шаг 1. Вычислить стадию k₁ по формуле (2).

Шаг 2. Вычислить стадию k₂ по формуле (2).

Шаг 3. Вычислить оценку локальной ошибки δ_n,1 по первой формуле (23).

Шаг 4. Вычислить значение q по формуле

                                               .                                                                           (34)

Шаг 5. Если  q<1, то h положить q·h и перейти на шаг 2 (возврат).

Шаг 6. Вычислить стадии k_i, 3 ≤ i ≤ m по формулам (2).

Шаг 7. Вычислить значение f(y_n+1).

Шаг 8. Вычислить оценку локальной ошибки δ_n,1 по второй формуле (23).

Шаг 9. Перевычислить значение q по формуле (34).

Шаг 10. Если  q<1, то h присвоить q·h и перейти на шаг 2 (возврат).

Шаг 11. Вычислить приближенное решение y_n+1 по схеме (2).

Шаг 12. Вычислить значение v_n по формуле (32).

Шаг 13. Вычислить значение r по формуле r×v_n=γ_m,1.

Шаг 14. Вычислить прогнозируемый шаг h_n+1 по формуле

                                   .

Шаг 15. Если q×v_n>γ_m,1 и m<M, то положить m равным (m+1).

Шаг 16. Если m>3 и q×v_n<γ_m–1,1, то положить m равным (m–1).

Заключение

С помощью описанной схемы получены коэффициенты многочленов устойчивости до степени m=27, и построены соответствующие методы первого порядка точности. Известно, что максимальная длина интервала устойчивости m-стадийного метода типа Рунге-Кутта первого порядка точности равна 2m². В результате на каждое вычисление правой части приходится 2m единиц длины интервала устойчивости. Это означает, что если шаг ограничен по устойчивости, то с ростом m эффективность метода возрастает. Построенный здесь алгоритм интегрирования жестких задач с переменным числом стадий на основе методов первого порядка позволяет существенно повысить эффективность расчетов на участке установления решения, где шаг ограничен по устойчивости.

Работа выполнена при финансовой поддержке РФФИ (проект 14-01-00047).

 

Литература:

1.      Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи: монография / Э. Хайрер, Г. Ваннер. – М.: Мир, 1999. – 685 c.

2.      Новиков Е.А. Явные методы для жестких систем: монография / Е.А.Новиков. – Новосибирск: Наука, 1997. – 197 с.

3.      Новиков A.Е. L-устойчивый (2,1)-метод решения жестких неавтономных задач / А.Е. Новиков, Е.А.Новиков // Вычислительные технологии. – 2008. – №13. – С. 477−482.

4.      Новиков Е.А. Алгоритм интегрирования жестких систем на основе (m,k)-метода второго порядка точности с численным вычислением матрицы Якоби / Е.А. Новиков, Ю.А. Шитов // Красноярск: Препринт ВЦ СО АН СССР №20. – 1988. – 23 c.

5.      Новиков Е.А. Численный алгоритм построения многочленов устойчивости методов первого порядка / Е.А.Новиков, М.В. Рыбков // Вестник Бурятского государственного университета. – 2014. – № 9-2. – С. 80–85.