Научный журнал
Фундаментальные исследования
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,749

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

Новиков Е.А. 1 Ващенко Г.В. 2
1 ФГБУН «Институт вычислительного моделирования СО РАН»
2 ФГБОУ ВПО «Сибирский государственный технологический университет»
Представлен L-устойчивый метод третьего порядка точности для численного решения жестких задач. Этот метод достаточно прост в реализации и обладает хорошими свойствами по точности. Эффективность предлагаемого метода может быть повышена за счет использования алгоритмов с контролем точности вычислений на шаге. Для реализации пошагового контроля нами построен критерий, в основе которого лежит оценка аналога глобальной ошибки. Эта оценка осуществляется с привлечением ранее вычисленных стадий численной схемы, что обеспечивает возможность выбора величины шага интегрирования фактически без увеличения вычислительных затрат. Сформулированы и представлены последовательный алгоритм и параллельная версия MPI-алгоритма с контролем точности вычислений на шаге. Данный метод входит в семейство (m, k)-методов, введенных в работе [5]. В (m,k)-методах стадия метода не связывается с обязательным вычислением правой части исходной задачи. За счет этого их свойства могут быть улучшены. Данные схемы можно рассматривать как способ реализации неявных методов типа Рунге ‒ Кутты. Важно, что при такой реализации не требуются итерации метода Ньютона.
контроль точности
переменный шаг
(3
2)-метод
параллельный MPI-алгоритм
1. Ващенко Г.В. Построение соотношения изоэффективности параллельного (3,2)-метода для численного решения задачи Коши // Современное состояние естественных и технических наук: мат. X Международной конф. – М.: Изд-во «Спутник+», 2013. – C. 36–41.
2. Демидов Г.В., Юматова Л.А. Исследование некоторых аппроксимаций в связи с А–устойчивостью полуявных методов // Численные методы механики сплошной среды. – 1977. – Т. 8. – № 3. – С. 68–79.
3. Новиков В.А., Новиков Е.А., Юматова Л.А. Замораживание матрицы Якоби в методе типа Розенброка второго порядка точности // ЖВМ и МФ. – 1987. – Т. 27. – № 3. – С. 385–390.
4. Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем – Новосибирск: Изд-во НГТУ, 2012. – 450 с.
5. Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые безытерационные методы решения жестких систем // ДАН СССР. – 1988. – Т. 301. – № 6. – С. 1310–1314.
6. Демидов Г.В., Юматова Л.А. Исследование некоторых аппроксимаций в связи с А–устойчивостью полуявных методов // Численные методы механики сплошной среды. – 1977. – Т.8. – № 3. – С. 68–79.
7. Новиков Е.А. Явные методы для жестких систем. – Новосибирск: Наука, 1997. – 197 с.
8. Новиков Е.А. Ващенко Г.В. Design of Parallel Algorithm (3,2)-method for Stiff Problems //Second International Conference on Cluster Computing (CC 2013), Lviv, Ukraine, June 3–5, Proceedings, p.165, 2013.
9. Demmel J. Applied numerical linear algebra. SIAM, 1997.
10. Hairer E. and Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems. Springer-Verlag, second revised edition, 1996.
11. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. – 1963. – № 5. – Р. 329–330.
12. Shampine L.M. Implementation of Rosenbrock methods // ACM Transaction on Mathematical Software. – 1982. – Vol. 8, № 5. – Р. 93–113.

Для численного решения жестких систем обыкновенных дифференциальных уравнений обычно применяются L-устойчивые методы [9]. При реализации таких численных схем на каждом шаге несколько раз решается линейная система алгебраических уравнений с применением LU-разложения матрицы Якоби. При большой размерности исходной задачи общие вычислительные затраты фактически полностью определяются временем декомпозиции данной матрицы. Сокращения затрат достигают за счет применения одной матрицы на нескольких шагах [3, 4, 9]. Наиболее естественно это осуществляется в итерационных методах, где данная матрица только определяет скорость сходимости итерационного процесса [9]. Для безытерационных методов типа Розенброка [10] и их модификаций [5, 9] вопрос о замораживании матрицы Якоби более сложный. В таких методах эта матрица включена в численную схему, и поэтому ее аппроксимация приводит к потере порядка точности численной формулы. Максимальный порядок методов типа Розенброка с замораживанием матрицы Якоби равен двум [3]. Безытерационные методы просты с точки зрения реализации и, как следствие, привлекательны для вычислителей. Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и L-устойчивых методов с автоматическим выбором численной схемы. В этом случае эффективность алгоритма может быть повышена за счет расчета переходного участка, соответствующего максимальному собственному числу матрицы Якоби, явным методом. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [6, 11]. Следует отметить, что применение таких комбинированных алгоритмов полностью не снимает проблему замораживания матрицы Якоби, потому что явным методом можно просчитать, вообще говоря, погранслойное решение, соответствующее максимальному собственному числу матрицы Якоби. В [5] предложен новый класс одношаговых численных схем, которые были названы (m, k)-методами. Они столь же просты при реализации, как и методы типа Розенброка, но обладают более хорошими свойствами точности и устойчивости. Более существенно, они достаточно просто реализуются с замораживанием матрицы Якоби.

Здесь разработан L-устойчивый (3,2)-метод третьего порядка точности для решения жестких задач. Построено неравенство для контроля точности вычислений, основанное на оценке аналога глобальной ошибки. Оценка осуществляется с привлечением ранее вычисленных стадий, что позволяет выбирать величину шага интегрирования фактически без увеличения вычислительных затрат. Сформулирован последовательный алгоритм и его параллельный аналог – MPI-алгоритм.

1. Класс (m, k)-методов решения жестких задач

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

y′ = f(y); y(t0) = y0; t0 ≤ t ≤ tk, (1)

где y и f – вещественные N-мерные вектор-функции; t – независимая переменная. Пусть Z – множество целых чисел, и заданы числа m, k ∈ Z, k ≤ m. Обозначим через Mm множество чисел {i ∈ Z, 1 ≤ i ≤ m}, а через Mk и Ji, 1 < i ≤ m подмножества из Mm вида

novikov001.wmf

novikov002.wmf

где 1 ≤ i ≤ m. Тогда класс (m, k)-методов записывается следующим образом [5]:

novikov003.wmf novikov004.wmf

novikov005.wmf i ∈ Mk, (2)

novikov006.wmf i ∈ Mm\Mk,

где ki, 1 ≤ im, - стадии метода; a, pi, βij, αij и cij - постоянные коэффициенты; h - шаг интегрирования; E - единичная матрица; fn = ∂f(yn)/∂y - матрица Якоби системы (1); - количество вычислений функции f на шаге; m - число стадий. На каждом шаге интегрирования осуществляются одно вычисление матрицы Якоби и одна декомпозиция матрицы Dn. Так как k и m полностью определяют затраты на шаг, а набор чисел m1, ..., mk из множества Mk только распределяет их внутри шага, то методы типа (2) названы (m, k)-методами. Основное отличие приведенных методов (2) от существующих численных формул состоит в том, что в них стадия метода не связывается с обязательным вычислением правой части исходной задачи (1), за счет этого свойства методов могут быть улучшены.

При k = m и αij = cij = 0 схемы (2) совпадают с методами типа Розенброка [10], а при k = m и αij = 0 – с ROW-методами [9]. В отличие от ROW-методов в численных формулах (2) более точно определены затраты на шаг интегрирования и более правильно описана область определения коэффициентов численных формул, что упрощает их исследование и делает их предпочтительнее. При рассмотрении методов такого типа в основном изучался случай k = m, то есть когда число стадий и количество вычислений правой части задачи (1) совпадают. В этом случае k-стадийную схему (2) можно поставить в соответствие k-стадийной полуявной формуле типа Рунге ‒ Кутты, при реализации которой на каждом шаге используется одна матрица размерности N. Относительно таких численных формул известно, что нельзя построить k-стадийную схему выше (k + 1)-го порядка точности, причем схема максимального порядка A-устойчивая. Если рассматривать схемы (2) при m > k, то можно построить численные схемы более высокого порядка точности [5].

2. Исследование (3, 2)-метода

Рассмотрим численную формулу вида

novikov007.wmf (3)

Разлагая стадии ki, 1 ≤ i ≤ 3 в ряды Тейлора и подставляя в первую формулу (3), получим ряд Тейлора для приближенного решения yn + 1. Полагая yn = y(tn) и сравнивая ряды для точного и приближенного решений, получим условия третьего порядка точности схемы (3), то есть

novikov008.wmf novikov009.wmf

novikov010.wmf novikov011.wmf

Положим a, β31 и β32 свободными и исследуем эту систему на совместность. Получим

novikov012.wmf

novikov013.wmf (4)

novikov014.wmf novikov015.wmf

где β = β31 + β32. Исследуем устойчивость схемы (3) на линейном скалярном уравнении y′ = λy, где смысл комплексного числа λ, Re(λ) < 0, – некоторое собственное число матрицы Якоби задачи (1). Применяя (3) для решения этого уравнения, получим yn + 1 = Q(z)yn, где z = hλ, а функция устойчивости Q(z) записывается следующим образом:

novikov016.wmf

Из вида Q(z) следует, что для L-устойчивости схемы (3) необходимо выполнение соотношения

a2 - a(p1 + p3) + β31p3 = 0.

Подставляя сюда коэффициенты (4), получим уравнение

a3-3a2 + 1,5a - 1/6 = 0.

Далее, сравнивая представление приближенного и точного решений до членов с h4 включительно, видим, что слагаемые с элементарными дифференциалами f′′′f3 и f′′ff2 в главном члене локальной ошибки будут отсутствовать, если

novikov017.wmf

novikov018.wmf

Теперь отсюда и (4) получим набор коэффициентов

novikov019.wmf

novikov020.wmf

novikov021.wmf

novikov022.wmf

novikov023.wmf

novikov024.wmf,

при которых локальная ошибка δn,3 схемы (3) имеет вид

novikov025.wmf

где значение a определяется из условия L-устойчивости

a3 – 3a2 + 1,5a – 1/6 = 0.

Это уравнение имеет три вещественных корня. Согласно [2] требование A-устойчивости схемы (3) имеет вид 1/3 ≤ a ≤ 1,0685790, поэтому выбираем корень a = 0,435866521508459.

В жестких задачах поведение ошибки определяется элементарным дифференциалом f′3f [5]. Поэтому при построении оценки аналога глобальной ошибки будем учитывать только первое слагаемое в локальной ошибке. Для контроля точности вычислений и автоматического выбора величины шага интегрирования используем идею вложенных методов. Для этого рассмотрим двухстадийный метод (2) вида

novikov026.wmf

novikov027.wmf novikov028.wmf (5)

где yn вычислено по формуле (3). В численной формуле (5) применяются стадии метода (3), и поэтому она практически не приводит к увеличению вычислительных затрат. Нетрудно видеть, что при коэффициентах b1 = 0,5(4a – 1)/a и b2 = 0,5(1 – 2a)/a схема (5) имеет второй порядок точности, а ее локальная ошибка имеет вид

δ n,2 = (6a2 – 6a + 1)h3 f′2f/6 + O(h4).

Тогда в неравенстве для контроля точности можно применять оценку ошибки εn(jn) вида [5]

novikov029.wmf (6)

где

c = 4·|6a2 – 6a + 1|/|1 – 12a + 36a2 – 24a3| ≈ 3.

При jn = 1 оценка εn(jn) будет A-устойчивой, а при jn = 2 – L-устойчивой. Теперь неравенство для контроля точности имеет вид

novikov030.wmf 1 ≤ jn ≤ 2, (7)

где ε – требуемая точность интегрирования, а параметр jn выбирается с наименьшим значением, при котором выполняется неравенство (7). Норма ||ξ|| в (7) вычисляется по формуле

novikov031.wmf

В случае выполнения неравенства novikov032.wmf по i-й компоненте решения контролируется абсолютная ошибка vε, в противном случае контролируется относительная ошибка ε.

3. Последовательный алгоритм интегрирования и его параллельная версия

Запишем схему (3) в покомпонентной форме, имеем

novikov033.wmf (8)

где 1 ≤ i ≤ N, novikov034.wmf novikov035.wmf и novikov036.wmf есть элементы векторов приращений novikov037.wmf novikov038.wmf и novikov039.wmf novikov040.wmf novikov041.wmfnovikov042.wmf и novikov043.wmf есть элементы матрицы Dn и векторов g(n) , σ(n) , w(n) ;

novikov044.wmf

novikov045.wmf

novikov046.wmf

novikov047.wmf

при i = j и novikov048.wmf при i ≠ j. Здесь novikov049.wmf – элементы матрицы Якоби novikov050.wmf задачи (1), вычисленные на решении y(n), novikov051.wmf – элементы вектора правой части f(y). Применение LU-разложения [8] приводит к системам уравнений для нахождения стадий novikov052.wmf novikov053.wmf и novikov054.wmf, то есть

novikov055.wmf novikov056.wmf 1 ≤ i ≤ N;

novikov057.wmf novikov058.wmf 1 ≤ i ≤ N; (9)

novikov059.wmf novikov060.wmf 1 ≤ i ≤ N,

где novikov061.wmf при i ≤ j и novikov062.wmf при i > j.

Используя обозначения элементов обратной матрицы novikov063.wmf через novikov064.wmf, нормы ошибок вычисляем по формулам

novikov065.wmf (10)

novikov066.wmf (11)

Пусть приближение yn к решению y(tn) задачи (1) вычислено в точке tn с шагом hn. Тогда учитывая, что имеет место novikov067.wmf 1 ≤ jn ≤ 2, последовательный алгоритм интегрирования формулируется следующим образом.

Шаг 1, 2. Вычислить матрицу Якоби novikov068.wmf и сформировать матрицу novikov069.wmf

Шаг 3. Выполнить декомпозицию матрицы novikov070.wmf

Шаг 4, 5. Вычислить

novikov071.wmf

novikov072.wmf

Шаг 6, 7. Вычислить

novikov073.wmf

novikov074.wmf

Шаг 8. Вычислить

novikov075.wmf

Шаг 9, 10. Вычислить норму ошибки novikov076.wmf по формуле (10) и q1 по формуле

novikov077.wmf.

Шаг 11. Если q1 < 1, то есть требуемая точность не достигнута, то вычисляется novikov078.wmf по формуле (11). В противном случае novikov079.wmf полагается равным novikov080.wmf.

Шаг 12. Вычислить значение параметра q2 по формуле novikov081.wmf.

Шаг 13. Если q2 < 1, то hn полагается равным q2hn и возврат на шаг 2.

Шаг 14. Вычислить приближение к решению в точке tn + 1 по формуле (8), то есть

novikov082.wmf

Шаг 15. Вычислить значение hn + 1 по формуле hn + 1 = min(q1, q2)hn.

Шаг 16. Выполнить следующий шаг интегрирования.

Замечание

При численном вычислении матрицы Якоби шаг численного дифференцирования ri, 1 ≤ i ≤ N, выбирается по формуле [6]

novikov083.wmf 1 ≤ j ≤ N,

где rmin – минимальный шаг численного дифференцирования, зависит от разрядности вычислительной системы. При расчетах с двойной точностью величину rmin следует принять равной 10–14. Тогда j-й столбец novikov084.wmf численной матрицы вычисляется по формуле

novikov085.wmf

то есть для задания матрицы требуется N вычислений правой части системы дифференциальных уравнений (1).

Если рассматривать алгоритм (3, 2)-метода (3), (8) как объект для распараллеливания, то его последовательный вариант выглядит как процесс вычисления векторов приращений novikov086.wmf, 1 ≤ i ≤ 3. При этом на каждом n-м шаге вычисления имеют последовательный порядок, novikov087.wmf При построении параллельного алгоритма необходимо сохранить этот порядок вычисления. Элементы каждого из векторов приращений получаются из решения систем линейных уравнений с одинаковой матрицей Dn и разными правыми частями.

Предположим, что размерность N системы (1) связана с размером p вычислительной системы соотношением N = s·p. Для Dn выберем блочно-строчную схему хранения. Тогда параллельный аналог (3,2)-метода (3) запишется в виде [7]

novikov088.wmf (12)

где 1 ≤ i ≤ p, novikov089.wmf.

Теперь сформулируем параллельный алгоритм вычисления приближенного решения с контролем точности вычислений. Для контроля точности в (12) используем процессор pr(1). Пусть известно приближение y(n) в точке tn с шагом hn. Тогда для вычисления y(n + 1) в точке t n + 1 справедлив параллельный алгоритм, в котором на каждом процессоре pr(j) формируется своя sj-я часть матрицы Dn, векторов novikov090.wmf novikov091.wmf novikov092.wmf и вектора решения novikov093.wmf.

Шаг 1. В каждом pr(j), 1 ≤ i ≤ p, novikov094.wmf:

а) выполнить recv novikov095.wmf

б) вычислить novikov096.wmf

в) вычислить матрицу Якоби.

Шаг 2, 3. Сформировать матрицу novikov097.wmf и разложить

novikov098.wmf.

Шаг 4, 5. Вычислить

novikov099.wmf

novikov100.wmf.

Шаг 6. В каждом pr(j), 1 ≤ i ≤ p, novikov101.wmf:

а) выполнить

recv novikov102.wmf

б) вычислить

novikov103.wmf

в) выполнить send novikov104.wmf.

Шаг 7. В каждом pr(j), 1 ≤ j ≤ p, novikov105.wmf:

а) выполнить recv novikov106.wmf

б) сформировать novikov107.wmf

в) вычислить novikov108.wmf и novikov109.wmf

Шаг 8. Вычислить

novikov110.wmf.

Шаг 9. В каждом pr(j), 1 ≤ j ≤ p, novikov111.wmf:

а) вычислить

novikov112.wmf

б) вычислить

novikov113.wmf

в) вычислить

novikov114.wmf

г) вычислить

novikov115.wmf

д) send novikov116.wmf

Шаг 10. В pr(1):

а) выполнить

recv novikov117.wmf

выполняется последовательность действий по контролю точности и вычисление значения hn;

б) выполнить send (hn, rp; 1, …, p). При rp = 1 – возврат на шаг 2, при rp = 0 – переход на шаг 11.

Шаг 11. В каждом pr(j), 1 ≤ j ≤ p, novikov118.wmf:

а) выполнить recv (hn, rp; 1);

б) вычислить novikov119.wmf

в) выполнить send novikov120.wmf.

Шаг 12. В pr(1):

а) вычислить h n + 1 по формуле

h n + 1 = min(q1, q2)hn;

б) send (h n + 1; 1, …, p).

Шаг 13. Выполнить следующий шаг интегрирования.

Замечание

Представленный алгоритм является параллельно-последовательным. В нем не учтены фрагменты, относящиеся к вычислению правой части (1) и матрицы Якоби, а также обращение матрицы Якоби, необходимое при невыполнении неравенства по точности. Для LU-разложения используется частичный выбор по столбцу.

В (m, k)-методах стадия метода не связывается с обязательным вычислением правой части исходной задачи. За счет этого их свойства могут быть улучшены. Данные схемы можно рассматривать как способ реализации неявных методов типа Рунге – Кутты. Важно, что при такой реализации не требуются итерации метода Ньютона, а все проблемы решаются за счет выбора шага интегрирования. Построена параллельная MPI-версия алгоритма интегрирования переменного шага, ориентированная на кластерные системы и топологию полного графа. В [1] построено соотношение изоэффективности, которое может быть использовано для сравнения различных параллельных алгоритмов решения одной и той же задачи Коши на основе (3,2)-метода, а также подходов к выбору и построению алгоритмов вычисления матрицы Якоби и алгоритмов ее факторизации. В особенности это относится к оценке коммуникационных затрат при организации межпроцессорных обменов.

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

Рецензенты:

Белолипецкий В.М., д.ф.-м.н., профессор, главный научный сотрудник, ФГБУН «Институт вычислительного моделирования» СО РАН, г. Красноярск;

Плотников С.М., д.т.н., профессор, ФГБОУ ВПО «Сибирский государственный технологический университет», г. Красноярск.

Работа поступила в редакцию 10.07.2014.


Библиографическая ссылка

Новиков Е.А., Ващенко Г.В. L-УСТОЙЧИВЫЙ МЕТОД ТРЕТЬЕГО ПОРЯДКА ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ ЖЕСТКИХ ЗАДАЧ // Фундаментальные исследования. – 2014. – № 9-4. – С. 734-740;
URL: https://fundamental-research.ru/ru/article/view?id=34917 (дата обращения: 29.11.2021).

Предлагаем вашему вниманию журналы, издающиеся в издательстве «Академия Естествознания»
(Высокий импакт-фактор РИНЦ, тематика журналов охватывает все научные направления)

«Фундаментальные исследования» список ВАК ИФ РИНЦ = 1.074