Software & Systems

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

Известные абсолютно устойчивые схемы фактически аппроксимируют только гиперболическое уравнение, а параболическое уравнение аппроксимируют при t/h ® 0, то есть при выполнении почти столь же жесткого условия, как и для классической явной схемы В.К. Саульева и Дюффорта–Франкела [1, 2]. Существуют абсолютно устойчивые схемы со счетом по явным формулам, называемые «классики» [3, 4]. В работе [3] доказывается сходимость метода с первым порядком в предположении t/h2 = const. Таким образом, абсолютно устойчивые схемы Саульева, Дюффорта–Франкела и «классики» аппроксимируют уравнение теплопроводности и сходятся при выполнении почти столь же жестких условий, как и для классической явной схемы [5–7]. Обобщить закон Фурье, добавив производную потока по времени, умноженную на малый параметр, предлагается в [8]. При этом уравнение теплопроводности становится сингулярно возмущенным гиперболическим уравнением, а условие устойчивости явной схемы для него менее жестким. В работе [9] приведены оценки разности решений задач для параболического и гиперболического уравнений с малым параметром. В [10] предложен алгоритм повышения устойчивости явных схем путем адаптивного сглаживания при помощи спектрально-сингулярного анализа. Однако в данном случае на каждом временном шаге необходимо находить сингулярное разложение траекторной матрицы решения, что может существенно замедлять вычисления. Поэтому при создании программ для многопроцессорных (многоядерных) систем желательно иметь явные схемы, аппроксимирующие уравнение теплопроводности и обладающие абсолютной устойчивостью. Счет по таким схемам не должен быть существенно медленнее счета по простым неявным схемам.

В данной работе для достижения безусловной устойчивости счета предлагается использовать регуляризацию решения, полученного на каждом временном шаге по классической явной разностной схеме, при помощи алгоритма Кули–Тьюки быстрого преобразования Фурье [11, 12]. При этом трудоемкость предлагаемых алгоритмов благодаря алгоритму Кули–Тьюки примерно равна трудоемкости метода прогонки.

Описание метода. Задача Дирихле

Рассмотрим начально-краевую задачу для простейшего уравнения теплопроводности:

 x Î [0, 1], u(0, t) = u(1, t) = 0,

u(x, 0) = f(x).                                           (1)

Ее аналитическое решение:

,

где .

Введем равномерную пространственную сетку  с шагом h = 1/N, содержащую N + 1 узел.

Аппроксируем уравнение (1) на этой сетке в виде

, ,

u0 = uN = 0, ui(0) = f(xi),

и запишем решение этой задачи:

xi = ih, i Î 0 ¸ N,

, .

Пусть , тогда  и , то есть только первые собственные элементы разностной аппроксимации аппроксимируют дифференциальный оператор и поэтому необходимо правильно рассчитывать эволюцию во времени только нескольких первых низкочастотных компонент и подавлять рост высокочастотных. Это обстоятельство, вероятно, впервые отмечено в [13]. Решением задачи (1) является ряд Фурье, и, согласно А.Н. Тихонову [14, 15], для получения устойчивого решения задачи суммирования ряда Фурье следует заменить исходный ряд на регуляризированный:

,

где ; a – параметр регуляризации. Легко заметить, что устойчивое суммирование ряда Фурье по Тихонову сводится к подавлению высших гармоник с ростом их номера.

Рассмотрим с этой точки зрения две классические разностные схемы для задачи (1) первого порядка аппроксимации по t и h:

условно устойчивую явную схему

                   (2)

и абсолютно устойчивую неявную схему

,               (3)

а также абсолютно устойчивую схему Кранка–Николсон (схему с весами) второго порядка аппроксимации по t и h:

           (4)

В момент времени t = nt полагаем

.

Тогда  и для схем (2)–(4) выполняются следующие соотношения:

,                    (5)

,                          (6)

.                      (7)

Из формулы (5) следует условие устойчивости явной схемы

Построим три явные схемы, отличающиеся алгоритмами регуляризации. Предварительно введем обозначения для дискретных прямого и обратного регуляризованного преобразований Фурье:

,          (8)

.                (9)

Схема 1. Введем параметр  определяющий качество аппроксимации собственного значения дифференциальной задачи собственным значением разностной.

Поскольку  для устойчивости гармоник с k £ k0 необходимо выбрать временной шаг интегрирования , где k0 определяется из неравенства  то есть

.                                     (10)

Квадратные скобки в формуле (10) означают целую часть числа.

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

.

Затем находим дискретное преобразование Фурье:  и, используя регуляризующие множители

вычисляем решение на n+1-м временном шаге: .

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

Схема 2. На каждом временном шаге строим дискретное преобразование Фурье разностной аппроксимации второй производной:

, ,

вычисляем регуляризованную вторую производную  с регуляризующими множителями  и в итоге находим решение на n+1-м временном слое:

.                                        (11)

Данный алгоритм является явным, эквивалентным неявной схеме и абсолютно устойчивым. Действительно, если положить

 и ,                  (12)

то при непосредственной подстановке (12) в (11) легко убедиться, что результатом алгоритма является соотношение (6).

Схема 3. На каждом временном слое сначала находим промежуточное решение:

, а затем преобразование Фурье промежуточного решения:  После этого с помощью обратного регуляризованного преобразования Фурье вычисляем решение  на n+1-м временном слое. В данном случае регуляризирующие множители определяются по формуле

.

Нетрудно убедиться в том, что при  и  из алгоритма следует соотношение , совпадающее с соотношением (7).

Алгоритм явной схемы является явным, абсолютно устойчивым и, хотя промежуточное решение не аппроксимирует исходное уравнение, эквивалентным схеме Кранка–Николсон.

Пример расчета. Задача Дирихле. Рассмотрим модельный пример задачи Дирихле:

,

u(x, 0) = 0, u(0, t) = u(l, t) = 0.

Задача имеет аналитическое решение:

.

На рисунке 1а приведены графики аналитического u(x) и численного ui решений (основ-ная вертикальная ось) и график относительной ошибки d = … в % (вспомогательная ось). Для расчета ui используется схема 3, временной шаг интегрирования в 105 раз превосходит максимальный шаг, определяемый из классического условия устойчивости явной схемы, число пространственных узлов N = 256, временных шагов 130.

В данном случае аналитическое и численное решения практически совпадают.

На рисунке 1б приведены графики аналитического и численного решений (основная вертикальная ось) и график относительной ошибки в % (вспомогательная ось). Для счета использована схема 3, временной шаг интегрирования в 107 раз превосходит максимальный шаг, вытекающий из классического условия устойчивости явной схемы. В этом примере число пространственных узлов N = 256, временных шагов 130. Графики аналитического и численного решений резко отличаются, относительная ошибка везде составляет порядка 50 %.

На рисунке 1 представлены результаты численно устойчивого счета, однако разностная схема с временным шагом в 107*h2/2 плохо аппроксимирует дифференциальную задачу. В таблицах 1–3 приведены некоторые результаты сравнения численного и аналитического решений для схем 1–3. Для сравнения найдем максимальную относительную ошибку на сетке в момент T (в %):

       (13)

здесь u(xi, T) – аналитическое, а  – численное решение в момент времени Т в точке xi. В таблице 1 приведены результаты счета по схеме 1. Здесь l = 2, N – размерность пространственной сетки, t/t0 – отношение шага интегрирования по времени к максимальному временному шагу, допускаемому условием устойчивости, nstep – количество выполненных временных шагов. В таблице 2 отражены некоторые результаты счета по схеме 2 при N = 256, в таблице 3 – результаты счета по схеме 3 при N = 256, nstep = 130.

Приведенные в таблицах примеры демонстрируют абсолютную устойчивость явных разностных схем 2 и 3 и значительное ослабление условия устойчивости для схемы 1.

Модификация метода для задачи Неймана. Рассмотрим решение краевой задачи с однородными условиями Неймана на границах х = 0 и х = 1. В этом случае можно использовать построенные выше схемы, в которых достаточно заменить синус-преобразование Фурье на косинус-преобразование Фурье с теми же регуляризирующими множителями. Несколько более сложная ситуация возникает, если на одном конце, например, при х = 0, задано условие Дирихле, а при x = l – условие Неймана. Вместо задачи (1) рассмотрим задачу , x Î [0, 1], u(0, t) = 0,  u(x, 0) = f(x).

Ее аналитическое решение:

Схемы 1–3 можно применить, если использовать дискретное преобразование Фурье по функциям  При этом, однако, неудобно применять алгоритм Кули–Тьюки быстрого преобразования Фурье. Гораздо удобнее симметричным образом продолжить решение и свести задачу к задаче Дирихле на отрезке удвоенной длины. В этом случае изменения в схеме 1 коснутся дискретного преобразования Фурье (8), (9), которое примет вид

,

.

Для нахождения решения на каждом n+1-м временном слое сначала вычисляем промежуточное решение на сетке из 2N узлов:

,

,

Далее при помощи алгоритма быстрого преобразования Фурье находим дискретное преобразование Фурье  и, используя регуляризирующие множители

 вычисляем , где i принимает значение от 1 до N. Значение k0 определяется формулой (10). Граничное условие  при x = l выполняется на каждом временном слое, так как .

Схемы 2 и 3 модифицируются аналогичным образом.

Пример расчета. Рассмотрим модельный пример: , u(0, t) = 0,

, u(x, 0) = 0.

Эта задача имеет аналитическое решение:

.

На рисунке 2 приведены графики аналитического и численного решений этой задачи (основная вертикальная ось) и график относительной ошибки в % (вспомогательная ось). Для счета также использована схема 3, временной шаг интегрирования в 105 раз превосходит максимальный шаг, вытекающий из условия устойчивости явной схемы. В этом примере 256 пространственных узлов и 156 временных шагов. Графики аналитического и численного решений практически сливаются.

Некоторые результаты счета этой задачи по схеме 1 приведены в таблице 4, по схеме 2 – в таблице 5, по схеме 3 – в таблице 6. Число пространственных узлов N = 256, временных шагов 156.

Упрощение регуляризации. Данные, приведенные в таблицах 1–6, свидетельствуют об устойчивости счета по схемам 1–3. Заменой пространственных переменных реальная задача может быть сведена к приведенным схемам. Для многих квазилинейных уравнений эти замены необходимо делать на каждом вре-менном шаге и при этом пересчитывать регуляризирующие множители. Для упрощения процедуры пересчета воспользуемся тем, что  и экстраполяция для всех k соотношения , справедливого при , приведет к более сильному подавлению высших гармоник разложения Фурье. Поэтому вместо регуляризирующего множителя  полагаем .

Таким образом, используем регуляризирующий множитель Тихонова с параметром регуляризации , где l = h*N – длина интервала.

Пример расчета приведен в таблицах (http://www.swsys.ru/uploaded/image/2023-2/2023-2-dop/12.jpg), где сравниваются численное и аналитическое решения при использовании последних регуляризирующих множителей.

Заключение

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

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

Список литературы

1. Саульев В.К. Интегрирование уравнений параболического типа методом сеток. М.: Физматлит, 1960. 320 с.
2. Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.
3. Gourlay A.R. Hopscotch: A fast second order partial differential equation solver. IMA J. of Applied Math., 1970, vol. 6, no. 4, pp. 375–390.
4. Gourlay A.R., McGuire G.R. General hopscotch algorithm for the numerical solution of partial differential equation. IMA J. of Applied Math., 1971, vol. 7, no. 2, pp. 216–227.
5. Степанов М.М., Савельев С.К. Численные методы в ракетостроении. СПб, 2019. 211 с.
6. Мазо А.Б. Вычислительная гидродинамика. Ч. 1. Математические модели, сетки и сеточные схемы. Казань, 2018. 166 с.
7. Аристова Е.Н., Лобанов А.И. Практические занятия по вычислительной математике в МФТИ. Ч. 2. М.: МФТИ, 2015. 310 с.
8. Глотов В.Ю., Головизнин В.М., Четверушкин Б.Н. Балансно-характеристические разностные схемы для уравнений параболического типа // Матем. моделирование. 2020. Т. 32. № 4. С. 94–106. doi: 10.20948/mm-2020-04-07.
9. Репин С.И., Четверушкин Б.Н. Оценка разности приближенных решений задачи для параболического уравнения и гиперболического уравнения с малым параметром // Докл. Академии наук. 2013. Т. 451. № 3. С. 255–258. doi: 10.7868/S0869565213210056.
10. Бахмутский М.Л. Интегрирование нелинейных параболических уравнений при помощи явных разностных схем с адаптивным сглаживанием // Вестн. кибернетики. 2015. № 4. С. 72–82.
11. Cooley J.W., Tukey J.W. An algorithm for the machine calculation of complex Fourier series. Math. of Computation, 1965, vol. 19, no. 2, pp. 297–301.
12. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes. The Art of Scientific Computing. NY, Cambridge University Press, 2007, 1256 p.
13. Люстерник Л.А. О разностных аппроксимациях оператора Лапласа // УМН. 1954. Т. 9. № 2. С. 3–55.
14. Тихонов Н.А. Об устойчивых методах суммирования рядов Фурье // Докл. АН СССР. 1964. Т. 156. № 2. С. 268–271.
15. Маслаков М.Л., Егоров В.В. Регуляризация рядов Фурье с приближенными коэффициентами для задачи оценки плотности распределения вероятностей фаз // СибЖВМ. 2022. Т. 25. № 2. С. 157–171. doi: 10.15372/SJNM20220205.

Reference List

1. Saulev, V.K. (1960) Integration of Equations of Parabolic Type by the Method of Grids, Moscow, 320 p. (in Russ.).

2. Roach, P. (1980) Computational Fluid Dynamics, Moscow, 616 p. (in Russ.).

3. Gourlay, A.R. (1970) ‘Hopscotch: A fast second order partial differential equation solver’, IMA J. of Applied Math., 6(4), pp. 375–390.

4. Gourlay, A.R., McGuire, G.R. (1971) ‘General hopscotch algorithm for the numerical solution of partial differential equation’, IMA J. of Applied Math., 7(2), pp. 216–227.

5. Stepanov, M.M., Savelev, S.K. (2019) Numerical Methods for Rocket Building, St. Petersburg, 211 p. (in Russ.).

6. Mazo, A.B. (2018) Computational Fluid Dynamics. P. 1. Mathematical Models, Grids and Grid Schemes, Kazan, 166 p. (in Russ.).

7. Aristova, E.N., Lobanov, A.I. (2015) Practical Classes on Computational Math. in MPTI. P. 2, Moscow, 310 p. (in Russ.).

8. Glotov, V.Yu., Goloviznin, V.M., Chetverushkin, B.N. (2020) ‘Balance and characteristic finite difference schemes for equations of the parabolic type’, Math. Models Comput. Simul., 12(6), pp. 981–989. doi: 10.1134/S2070048220060095.

9. Repin, S.I., Chetverushkin, B.N. (2013) ‘Estimating the difference between approximate solutions of a problem for a parabolic equation and a hyperbolic equation with a small parameter’, Dokl. Math., 451(3), pp. 255–258. doi: 10.7868/S0869565213210056 (in Russ.).

10. Bakhmutsky, M.L. (2015) ‘Integration nonlinear parabolic equation by explicit finite-difference scheme with adaptive smoothing’, Bull. of Cybernetics, (4), p.72–82 (in Russ.).

11. Cooley, J.W., Tukey, J.W. (1965) ‘An algorithm for the machine calculation of complex Fourier series’, Math. of Computation, 19(2), pp. 297–301.

12. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. (2007) Numerical Recipes. The Art of Scientific Computing, NY: Cambridge University Press, 1256 p.

13. Lyusternik, L.A. (1956) ‘On difference approximations of the Laplace operator’, Russian Math. Surveys, 9(2), pp. 3–55 (in Russ.).

14. Tikhonov, A.N. (1964) ‘On stable summability methods for Fourier series’, Dokl. Akad. Nauk SSSR, 156(2), pp. 268–271 (in Russ.).

15. Maslakov, M.L., Egorov, V.V. (2022) ‘Regularization of Fourier series with approximate coefficients for the problem of phase probability density function estimation’, SJNM, 25(2), pp. 157–171. doi: 10.15372/SJNM20220205.