Численное решение краевых задач для уравнения Пуассона методом быстрого преобразования Фурье с использованием параллельных вычислений Текст научной статьи по специальности «Математика»

DOI: 10.25702/^^2307-5252.2018.9.5.165-182 УДК 533.95 + 519.63

О. В. Мингалев, М. Н. Мельник

ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ ДЛЯ УРАВНЕНИЯ ПУАССОНА МЕТОДОМ БЫСТРОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ С ИСПОЛЬЗОВАНИЕМ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ

Аннотация

Для численного решения краевых задач для уравнения Пуассона в 2-мерном прямоугольнике и 3-мерном параллепипеде на регулярной сетке с большим числом узлов методом быстрого преобразования Фурье (БФП) разработаны несколько новых приемов, которые позволяют эффективно использовать параллельные вычисления как на ядрах центрального процессора, так и на графических процессорах (GPU). Создан набор программ для случая периодических граничных условий, а также для случаев однородных граничных условий Дирихле и Неймана. Эти программы дают решение с 4-м порядком точности и свободны от ограничения на число шагов сетки по каждому измерению в исходном методе БФП. Программы имеют простой и удобный интерфейс, а также максимально возможный уровень параллельности, и позволяют достаточно быстро решать указанные задачи с числом узлов сетки порядка 109 и более.

Ключевые слова:

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

O. V. Mingalev, M. N. Melnik

NUMERICAL SOLUTION OF BOUNDARY PROBLEMS FOR THE EQUATIONS OF POISSON BY FAST FOURIER TRANSFORM USING PARALLEL CALCULATIONS

Abstract

For the numerical solution of boundary value problems for the Poisson equation in a 2D and a 3D on a regular grid with a large number of nodes, the Fast Fourier Transform (FFT) method has developed several new techniques that make it possible to efficiently use parallel computing on both the cores CPU and on graphics processing units (GPU). A set of programs has been created for the case of periodic boundary conditions, as well as for the cases of homogeneous Dirichlet and Neumann boundary conditions. These programs have a solution 4th order of accuracy and haven't the restriction on the number of grid steps for each dimension in the original FFT method. Programs have a simple and usable interface, as well as the highest possible level of parallelism, and allow you to quickly solve these problems with the number of grid nodes of the order of 109 or more.

Keywords:

Poisson's equation, of boundary value problems, the Fast Fourier Transform, Parallel calculation.

Введение

Проблема быстрого численного решения краевых задач для уравнения Пуассона (см., например, [1]) в 2-мерном прямоугольнике и 3-мерном

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

Для этих целей возможно использовать только такие численные методы, в которых основной объем вычислений может выполняться в параллельном режиме на графических процессорах (GPU). К числу таких относится метод быстрого преобразования Фурье (далее БФП) (см., например, [2-4]), который позволяет решать задачу с периодическими граничными условиями, а также задачи с однородными граничными условиями Дирихле и Неймана.

Отметим, что в платные библиотеки MKL и IMSL языка FORTRAN входят подпрограммы решения краевых задач для уравнения Пуассона в 2-мерном прямоугольнике и 3-мерном параллепипеде методом БФП, однако они имеют ряд недостатков, которые делают их непригодными для достижения указанных выше целей.

Метод быстрого преобразования Фурье допускает два варианта: со 2-м и 4-м порядком точности, и основан на алгоритме Кули-Туки (Кули-Тьюки) дискретного быстрого преобразования Фурье периодического 1 -мерного массива

с размерностью N = 2m , что дает соответствующее ограничение на числа Nxq ,

Nyо, Nzq шагов сетки по каждому измерению X, Y, Z вида

Nxo = 2mx, Nyo = 2my, Nzo = 2mz , (1)

где mx, my, mz e N - натуральные числа. Из этого ограничения вытекает условие

на соотношения размеров параллепипеда (прямоугольника) Lx , Ly , Lz по соответствующим измерениям X, Y, Z вида

Ly/Lx = 2My,х , Lz/Lx = 2Mz,x , My, e Z, M, e Z , (2)

где My, = my - mx, Mz, = mz - mx e Z - целые числа.

Это ограничение весьма неудобно для построения численных моделей крупномасштабных процессов в космической плазме. Отметим, что алгоритм Кули-Туки может быть обобщен также и для случая, когда размерность массива N имеет вид произведения степеней простых чисел 2, 3, 5, 7, то есть

дг _ 2>»2.3»% . 5Щ . i™! ^ где т2,т3,т5,т1 е Z+ - целые неотрицательные

числа. Однако в этом случае он существенно усложняется.

Для численного решения задачи Дирихле и задачи Неймана можно использовать два способа. Первый способ состоит в использовании соответственно дискретного

синус или косинус-преобразования Фурье. Второй способ состоит в сведении исходной задачи к периодической задаче в прямоугольной области с удвоенными размерами путем продолжения правой части специальным образом. Отметим, что быстрое дискретное синус и косинус-преобразование Фурье периодического 1 -мерного массива с размерностью N сводятся к обычному быстрому дискретному преобразованию Фурье периодического 1 -мерного массива с удвоенной размерностью 2N . При этом программы синус и косинус-преобразования существенно сложнее исходного преобразования Фурье. Кроме того, при таком сведении задачи Неймана к периодической автоматически достигается аппроксимация однородного граничного условия Неймана с тем же порядком точности, что и аппроксимация уравнения Пуассона. Поэтому второй способ для наших целей является более удобным, и в работе предложен прием сведения задачи Дирихле или Неймана с однородными граничными условиями в прямоугольнике [0;Ьх] х [0;Ьу ] или параллепипеде

[0; Ьх ] х [0; Ьу ] х [0; Ь2 ] к периодической задаче в аналогичной области с

удвоенными размерами по каждому измерению. При этом автоматически получается конечноразностная аппроксимация граничного условия Неймана того же порядка точности, что и аппроксимация уравнения Пуассона, в частности, аппроксимация 4-го порядка точности. Это является важным новым методическим результатом.

Также в работе предложен новый прием, который позволяет снять ограничение (1) на число шагов сетки по каждому измерению, то есть оставляет только условие (2), и состоит в следующем. Если условие (1) нарушено, то по исходной сетке с шагом й0 создается вспомогательная сетка с более мелким шагом И и числами шагов сетки по каждому измерению Nx, N у , N2 , которые удовлетворяют условию (1). При этом

числа тх, Шу, тг определяются через исходные числа шагов Nxo , N у 0 , N20 по формулам

Шх = N2 Nx0 ] +1> Шу =[^2 Nyo ] +1, т2 = [1о§2 Nzo ] +1, (3)

где через [а] обозначена целая часть действительного числа а. Для получения

дискретной задачи на вспомогательной сетке выполняем интерполяцию правой части уравнения с исходной сетки на вспомогательную сетку при помощи интерполяционного полинома Лагранжа с порядком на 2 большим используемого порядка аппроксимации уравнения. То есть в случае метода 2-го порядка точности используем полином Лагранжа 4-го порядка — аппроксимация по 4-м ближайшим узлам на измерение, а в случае метода 4-го порядка точности используем полином Лагранжа 6-го порядка — аппроксимация по 6-ти ближайшим узлам на измерение. После этого методом БФП решается задача на вспомогательной сетке, а потом полученное численное решение интерполируется обратно на исходную сетку полностью аналогично описанной интерполяции. Многочисленные тестовые расчеты показали высокую точность и эффективность этого приема.

В работе рассматриваются основные детали программы, которая выполняет следующие шаги.

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

2) Построение удовлетворяющей условию (1) вспомогательной сетки с меньшим шагом И , и интерполяция правой части уравнения на эту сетку, которую удобно выполнять в параллельном режиме.

3) Вычисление собственных чисел 1-мерной дискретной периодической задачи по каждому измерению на вспомогательной сетке.

4) Вычисление массива коэффициентов дискретного преобразования Фурье правой части на вспомогательной сетке. Эта первая из двух главных частей программы, которая содержит почти половину всех вычислений. Она выполняется в результате последовательного выполнения по каждому измерению 1-мерного быстрого преобразования Фурье для всех 1-мерных сечений массива правой части вдоль данного измерения. В ходе последнего процесса каждое 1 -мерное сечение массива правой части вдоль данного измерения независимо обрабатывается подпрограммой 1-мерного БФП. Поэтому эта часть программы удобна для выполнения в параллельном режиме как на ядрах центрального процессора, так и на графическом процессоре.

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

6) Вычисление сеточного массива решения при помощи обратного дискретного преобразования Фурье. Эта вторая из двух главных частей программы, которая по объему вычислений и методике выполнения аналогична прямому преобразованию Фурье, описанному в пункте 4).

7) Обратная интерполяция решения со вспомогательной сетки на исходную сетку, аналогичная прямой интерполяции в пункте 2).

8) Численное дифференцирование решения для нахождения соответствующего векторного поля, использующегося в модели.

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

1. Сведение задач Дирихле и Неймана к периодической задаче.

Обозначим в пространствах К2 и М3 векторы ортонормированного базиса декартовой системы координат через ех , еу и е2, а векторы в этих

пространствах будем рассматривать в виде:

2 ^ л- =хех + уеу еК , х=хех+ уеу+ .

Будем рассматривать 2-мерный открытый исходный прямоугольник П21 и

прямоугольник П2 2 с удвоенными размерами, которые определяются

формулами

П21 =(0;Ьх)х(0;Ьу) , П2,2 = {0,2ЬХ)х(0;2Ьу) (4)

а также аналогичные 3-мерные параллепипеды исходный П3 ^ и «удвоенный» П3 2 , которые определяются формулами

Пзд =(0;Ьх)х(0;Ьу)х(0;Ь2), =(0;2!х)х(0;2Ьу)х(0;2Ьг) . (5)

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

Рассмотрим краевые задачи Дирихле и Неймана с однородными граничными условиями и финитной правой частью для уравнения Пуассона

А и(х) = /(х) , х е Пи>1, (6)

где правая часть /(х) - заданная функция. Будем считать, что правая часть

У(х) е С^ПиД^ , то есть является достаточно гладкой (1 раз непрерывно

дифференцируема) на замкнутом прямоугольнике (параллепипеде) Пп,1 и финитной (равной нулю в некоторой окрестности границы ЗП^ ), где п = 2,3 -размерность пространства. Однородное граничное условие Дирихле имеет вид

и( х ) = 0, х е ЗПп1. (7)

Однородное граничное условие Неймана имеет вид

^ = 0, х е ЗП 1, (8)

он

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

часть /(х) с множества Пп,1 на удвоенный прямоугольник или параллепипед Пп,2 . Обозначим такое продолжение в случае задачи Дирихле через FD(х), а в случае задачи Неймана через FN( х) . Для продолжения FD(х) по каждой переменной используется нечетное продолжение относительно правого края Ьх , Ьу , Ьг исходного интервала. В 2-мерном случае продолжение FD(х) определяется формулами

^ ( х у ) = /( х, у) , ^ ( х + Ьх, уУ) = -/(Ьх- х, у ) , ^ ( ^ у + Ьу ) = -/( х, Ьу- у) , ^ ( х + Ьх, у + Ьу ) = /(Ьх- х, Ьу- у ),

(9)

( х, у ) е П2,1 = [0; Ьх ] х[0; Ьу ].

В 3-мерном случае продолжение FD( х) получается в результате последовательного применения следующих формул:

(х У, *) = у(х, У, *) , (х + Ьх, У, * ) = -/(Ьх- х, У, * ) Рв (x, У + Ьу, * ) = -/(х, Ьу- У, * ) ,

Рв (х + Ьх, У + ЬУ, *) =1 (Ьх - х, ЬУ - У, *)

при (х, У, * )е Пз,1 = [0; Ьх ] х [0; Ьу ] х [0; Ь* ], Рв ( x, У, * + Ь* ) = -/ ( х, У, Ь*- * )

(12)

при (х,У,*) е [0;2Ьх] х[0;2Ьу] х[0;Ь*].

Для продолжения Ру(х) по каждой переменной используется четное

продолжение относительно правого края исходного интервала. В 2-мерном случае оно о пределяется формулами

Рм (х у) = /(х, У) , Рм (х + Ьх, У) = у(Ьх- х, У) ,

Рм (х, У + Ьу ) = /(х, Ьу- у) , (х + Ьх, у + Ьу ) = /(Ьх- х, Ьу- у) , I (11)

(х, у) е П2,1 = [0; Ьх ] х [0; Ьу ], ^

а в 3-мерном случае получается в результате последовательного применения аналогичных формул:

Рм (х У, *) = У (х, У, *), Рм (х + Ьх, У, *) = У (Ьх- х, У, *), Рм ( x, У + Ьу, * ) = У ( х, Ьу- У, * ),

Рм (х + Ьх, У + Ьу, *) =1 (Ьх- х, Ьу- У, * ) при (х, у, * )е Пз,1 = [0; Ьх ] х [0; Ьу ] х [0; ^ ],

Рм ( х У, * + Ь* ) = у ( х, У, Ь*- * ) при ( х, у, * )е [0;2Ьх ] х [0;2Ьу ] х [0; Ь* ].

Рассмотрим в удвоенном прямоугольнике или параллепипеде Пи 2 периодические краевые задачи для уравнения Пуассона

ди (х ) = Р (х), х е Пи,2, (13)

с продолженными правыми частями Р( х) = х) и Р( х) = Рм( х) . Используя

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

Утверждение. Пусть /(х)еС^Пид) , где область Пи1 определяется формулами (4) и (5), а функции х) и Рм(х) являются продолжениями функции /(х) по формулам (9), (11) и (10), (12) соответственно на

определяемую формулами (4) и (5) удвоенную область Пи2 . Если ив (х) является классическим решением периодической краевой задачи (13) при

F (x) = Fd(x) , то его ограничение на Пид является классическим решением задачи Дирихле (6), (7). Если UN(x) является классическим решением периодической краевой задачи (13) при F (x) = FN(x) , то его ограничение на Пи,1 является классическим решением задачи Неймана (6), (8).

2. Интерполяция на вспомогательную сетку и обратно.

Рассмотрим вопрос построения вспомогательной сетки. Будем считать, что

размеры параллепипеда (прямоугольника) Lx , Ly , Lz удовлетворяют условию (2), и в задаче из физических условий введена сетка с шагом h = Lx/Nxq = LyjNyft = Lz/Nzq . Если числа шагов сетки по каждому

измерению Nx0 , Ny0 , Nz0 не удовлетворяют условию (1), то вводим

вспомогательную сетку с меньшим шагом h = Lx¡Nx = LyjNy = Lz/Nz , для

которой новые числа шагов сетки по каждому измерению Nx , N y , Nz

удовлетворяют условию (1). При этом числа mx, m.y, mz определяются через

исходные числа шагов Nx0 , Ny0, Nz0 по формулам (3). При помощи целых

векторов (мультииндексов) к = кх ex + ky ey и к = кх ex + ку ey при n = 2 , и

к = кх ex + ky ey + kz ez и к = kx ex + ky ey + kz ez при n = 3 узлы этих сеток обозначим как

r(k) = ho-k, f(k) = hk, к,к&Ъп, п = 2,3, (14)

Будем считать, что в исходной физической задаче задан массив значений Fh (к) = F(r (к)) правой части F(x) в узлах исходной сетки на замыкании

«удвоенной» области ПИ;2 . Для получения массива значений Fh (к^ = F^r (к^

правой части в узлах вспомогательной сетки выполняем интерполяцию с исходной сетки на вспомогательную при помощи интерполяционного полинома Лагранжа 6-го порядка — аппроксимация по 6-ти ближайшим узлам на измерение, то есть по 36 узлам при n = 2 и по 216 узлам при n = 3 . Применительно к рассматриваемому случаю эту формулу интерполяции при n = 2 удобно представить в форме

3 3

= Z Z Fhilhklh0] + Py

Px=-2Py=-2

^(К/М-ф-р! (15)

U S'x, Px + (px - 'x ) X

V x /

Í з /г.,- /, í .M. „ n

^y, Py +(py - ''y )

X

где через [Я ] и {Я} = Я - [Я] обозначены целая и дробная части действительного числа Я, через ■ к¡И{)^ = - кх /+ - ку /

\еу^гг и

р = рхех + ру еу е Z2 обозначены соответствующие векторы с целыми координатами, и через

П при I =р

= 1 п • ^ .

г [0 при гФ р

обозначен символ Кронекера. При п = 3, используя аналогичные обозначения, формулу интерполяции удобно представить в форме

3 3 3

¿и ¿^Ф^Ы^У

Рх =-2 Ру =-2 Рх =-2

Г з

п

V гх=-2

' 3

п

V гу=-2

( з

п

Vг*=-2

{{Ькх1К}-гх){1-81х^)

, рх +( Рх - гх )

({ЧЛЮМ ^

(16)

гу, Ру,

У, Ру

■(Ру -гУ)

, Р*

"( Р* - г* )

Обозначим через и^к} = и{г(к^ и Ок{к^ = и{г{к^ массивы значений решения соответственно в узлах исходной сетки и в узлах вспомогательной сетки,

Zn

Обратная интерполяция решения со вспомогательной сетки на исходную дается формулами, аналогичными (15) и (16):

3 3

Рх =-2 Ру =-2

3

п

Чгх =-2

({ИЛ/И}- гх )(1 , Рх )

Л

^ гх, Рх

(Рх - гХ )

(17)

^ ({И^у/И}-г'у)(1 у,Ру )

К =-2

гу, Ру

(Ру - гу )

п = 2.

х

3 3 3

U

■W=E Z ^Vb&dhok/hl+py,

px=-2 py=-2 px =-2

^ ({h0 kx/hb'x )(l -^x, px

V ix =-2

'x, px

( px - 'x )

П

iy=-2

3

П3

V 'z=-2

({ho ky/h}-iy )(

1 -5

'y > py

5iy , py +( py - ^ )

({ho k^h}-iz )(l -5'z, pz )'

5iz, pz +( pz - 'z )

, n = 3 .

(18)

3. Формулы для решения периодической задачи методом Фурье

Алгоритм решения дискретной краевой задачи методом БФП с логической точки зрения является полным аналогом решения соответствующей исходной «непрерывной» краевой задачи методом Фурье. Поэтому для отслеживания алгоритмической цепочки необходимо сначала рассмотреть формулы для решения периодической задачи (13) методом Фурье. В этих формулах используются собственные числа и функции 1 -мерной периодической задачи

d2V(x) _

dx

,2

= X-V (x), x е(0 ;2L ), V ( 0) = V (2L ), (19)

которые удобно обозначить следующим образом: \2

X[2í] = -(ot/Z) , F[2s](x) = cos(;rsx/Z), s = 0,1,2,... (s e Z+)

X[2s-1] = -fT 1 , V [2s-1] (x ) = sin ÍT

= 1,2,... (sgn)

(20)

и

Заметим, что собственные числа X[0] = 0 и X[4s - 3] = -(t(2s -1)/(2L)) являются однократными с собственными функциями соответственно У[0](х) = 1 и V[4s - 3](x) = sin (t(2s -1)x/( 2L)) , где s e N . Остальные собственные

числа являются двукратными: X[2s] = X[4s -1] = - (ts/L )2, и им соответствует пара собственных функций V[2s](x) = cos (rsx/L)

V[4s -1](x) = sin (Tisx/L) взаимно ортогональных в пространстве L2(0 ;2L) .

Собственные числа и функции для периодической задачи (13) удовлетворяют уравнению

AW (х ) = X -W( х), х е Пп2, (21)

и выражаются через определяемые формулами (20) собственные числа и функции 1-мерной периодической задачи (21) следующим образом. Обозначим

собственные числа и функции соответствующих 1-мерных задач вида (20) по переменным х , у и * как

2

М2р«] = -(лР«/ Ьо) ;

уЛ2Ра](х) = со*(лрах/1а), ра =0,1,2,... (ра&%+)

V« [2 Ра-1] ( х ) =

Х«[2 Ра-1] = -( Л Рах 1

V 2Ьа У

( \2 ЛРа

2Т

V 2Ьа У

ра= 1,2,... (^еН)

(22)

а = х, у, г.

Собственные числа и функции задачи (21) для п = 2 определяются формулами

х х г у "у

^1р](х) = Щ_рх,ру](х,у) = Ух1рх](х)-Уу1ру](у) ,

(23)

(24)

а для и = 3 определяются формулами

ЧР\ = К\Рх\ + \\Ру\+К\Рг\> Р = Рхех+Руеу+Ргег е

Обозначим через Р [р] коэффициенты Фурье правой части Р(х) уравнения Пуассона по определяемой формулами (22)-(24) ортогональной системе собственных функций {Ж[р](х)} задачи (21) в пространстве ¿2 (Пп2) , которые определяются формулой

| Р(х)Ж[р](х)^пх

Р[ Р] =

( Р ,Ж [ р])

Ь2(Пп,2 )

П

п,2

(Ж [ р],Ж [ р])

Ь2(Пп,2 )

| |Ж[р] (х)| ёпх

(25)

П

где через

(/1, у2 )Ь(Пи,2) = I у1(х)у2 (х)Л

П

п,2

обозначено скалярное произведение в пространстве Ь2 (Пп2) . Отметим, что условие разрешимости задачи (13) аналогично таковому для задачи Неймана и состоит в ортогональности правой части Р( х) первой собственной функции

Ж[0] (х) = 1 в пространстве Ь2 (Пп 2 ), то есть должно выполняться равенство:

>

2

+

| Е(х) йпх

= 0.

п

п,2

Алгоритм решения задачи (13) методом Фурье состоит из следующих шагов:

1) находим по формуле (25) коэффициенты Фурье Е [р] правой части по

определяемой формулами (22)-(24) системе собственных функций [ р]( х )| задачи (21) в пространстве ¿2 (Пп 2 );

2) находим коэффициенты Фурье и [ р] решения и (х) при

||р|| = Рх + Ру + Рг ^1 по формуле

и [р] = F [р]Д[р], (26)

и далее для удобства считаем, что и[0](х) = 0 ;

3) находим решение и (х) задачи (13) как сумму ряда Фурье

и ( х ) =

V V

и [Рх, Ру ] • Ух [Рх ] (х )-Уу [Ру ](у ) , п = 2.

/ , [Рх, Ру ] ' х [Рх ] Vх ) ' у [Ру ]

Рх =0 Ру =0

+ОТ +1» +ю

и ( х )=ЕХ^и [ Рх, Ру, Р* ]•

Рх =0 Ру =0 Рг =0 •Ух [Рх ]( х )• Уу [Ру ](у ) ' [Рг ](г ), п = 3.

(27)

4. Дискретная задача и ее решение методом Фурье

Алгоритм решения дискретной краевой задачи методом БФП полностью повторяет шаги изложенного выше алгоритма решения «непрерывной» краевой задачи (13), основанного на формулах (22)-(27), но при этом содержит ряд дополнительных технических действий.

Рассмотрим дискретный аналог задачи (19), то есть 1-мерную дискретную

периодическую задачу на собственные числа и векторы. Введем на отрезке [0 ;2Ь] сетку хк = к- И , к = 0,1,..,,2Ы , где к = Ь/Ы - шаг сетки. С учётом условия периодичности V (0) = V (2Ь) функция V (х) переходит в 2 N -мерный сеточный

вектор У/, е М2Л с координатами ¥ь(к) = ¥(хк), к = 0,1,.. .,2Ы-\.

В случае конечно-разностной аппроксимации 2-го порядка точности оператор 2-й производной заменяется разностным оператором 5 (см. [2-4]):

^Г- ¡И2^ (*) = (к + 1)~2п (*) + >г* (к- 0) • <28)

При этом 1-мерная периодическая задача (19) на собственные числа и функции

принимает вид задачи на собственные числа и векторы для оператора 5 (симметричной матрицы размера 2N х 2N):

8\(к) = Х-Г/г(к), к = 0,1,...,2^-1, (29)

где У^ (/с) определяется при всех А' е Ж из условия периодичности:

V (к ± 2Ы~) = Ук (к). В матричной форме задача (29) может быть записана как

А Ук = X ■ У11. где соответствующая оператору 5 матрица А размера 2 А' х 2 А' имеет вид:

A =

-21000 1 - 2 1 0 0 0 1 - 2 10

0 0 0 0 1 0 0 0 0 0 0 0 0 0 0

00000 00000 1 0 0 0 0

0 1-2 10 0 0 1-2 1 0 0 0 1-2

Нетрудно найти, что матрица A имеет (N +1) различных собственных чисел

X[s] = -2(1 + cos(xs/N)), s = 0,1,...,N, (30)

из которых два числа Xmin = Х[0] = —4 и Xmax = Х[N] = 0 однократны, а остальные (TV —l) чисел ).| v | при s = \,...,N — 1 являются двукратными. Учитывая, что формально при s = 0,1,.. .,2N верно равенство A[2iV — s] = Л[л],

собственные векторы £ М.2 V матрицы Á (и оператора б ) удобно

занумеровать так: s = 0,\,...,2N — \ .

При этом их координаты V[s](k) , k = 0,1,..., 2N — 1 определяются формулами

V[s](k) = cos(xsk/N) при s = 0,1,...,N , V [s](k) = sin {xsk/N) при s = N + 1,...,2N -1,

2 N-1

а для евклидовых норм собственных векторов |V[s]|2 = ^ |V[s](k)|2 верны

k=0

равенства:

V[0]|2 = V[N]|2 = N, V[s]f = N/2 при s = 1,..., N -1, N + 1,...,2N-1.

Рассмотрим теперь дискретный аналог задачи (13) в случае конечно-разностной аппроксимации 4-го порядка точности (см., например, [2,3]) на

(31)

$2 -2

определяемой формулой (14) вспомогательной сетке. Обозначим через ох, 5У и -2

б2 конечно-разностные аппроксимации 2-го порядка точности для операторов

д2 д2 д2

—т-, -т- и —- аналогичные (28):

дх ду дг

д2Ц(г(к)) 1 ;2 .

а = х, у, г.

В случае конечно-разностной аппроксимации с 4-м порядком точности (см., например, [2, 3]) оно при п = 2 переходит в дискретное уравнение:

ьЩк) = ±(иь(к + еа)-2иь(к) + иь(к-еа));

¥

г

2-2 1-2-21- ~ ~ ~ 1 /-2 -2

5*+5,, + -5*5,, |£/й(А) = ^(А) + —(5*+5,, |^(к), (32)

1 /-2 -2 \

а при п = 3 переходит в дискретное уравнение:

1

-2 -2 -2 1-2-2 1-2-2 1 -2 -2 ) ~ ~ бл+бд, +52 + -8*5>, + -5*52 + -8>,52 Шк(к) =

6 6 6 у

у (33)

~ ~ 1 /-2 -2 -2 \ ~ ~

Рассмотрим вспомогательную задачу на собственные числа и векторы стоящего в левой части уравнений (32) и (33) оператора. При п = 2 этот оператор действует в 4- Nx^ -мерном линейном пространстве матриц размера

(2^ х 2Nу ) с компонентами

ж(кх,ку), кх=0,\,...,2Мх-\, ку= ОД...,2^-1,

периодически продолженных по к на все Ъ2, которое обозначим как X2 . При п = 3 этот оператор действует в 8 • Nx• Ny • Nг -мерном линейном пространстве

3-мерных таблиц размера (2Nx х 2Nу х 2Nг) с компонентами

ж(кх,ку,кг), кх=0,\,...,2Мх-\, ку =0,\,...,2Ыу -1, К= 0,1,...,2^-1

также периодически продолженных по к на все 1?, которое обозначим через

. Уравнение для собственных чисел и векторов в зависимости от размерности задачи п имеет вид

ЬпЖ(к) = Х-Ж (к), п = 2,3,

- -2 -2 1 -2 -2 - -2 -2 -2 1 /-2 -2 -2 -2 -2 -2 Т>2 = 8^+8^+—8*8^ , Х)3 =8^+8^+82 + — (8x8^ +8*82 +8^8г 6 6 \

(34)

Собственные числа оператора можно представить через решение (30), (31) 1-мерной дискретной периодической задачи (29), которое дает собственные числа

и векторы операторов Ъх, 5 г и б2 в зависимости от размерности п в виде МР] = ^х [Рх ] + Ху [Ру ] + 1 Хх[Рх ] • Ху [Ру ] > п = 2,

Р] = Xх [Рх ] + Ху [Ру ] + Х* [Ру ] +

+ 1 (Х х [ Рх ] • Х у [ Ру ] + Х х [ Рх ] • Х * [ Ру ] + Х у [ Рх ] • Х * [ Ру ] ) . П = 3 =

(35)

где аналогично (30) и (22)

Х«[Рх ] = -2 (1 + СС8 (лРа/ N« ) ) , Р« = 0,1,...,2М«- 1,1 ^

а = х, у, *. \

При п = 2 взаимно ортогональные собственные векторы (матрицы) оператора образуют базис в указанном выше 4 • Мх • М^ -мерном линейном

пространстве периодических матриц и могут быть представлены в виде диады

$2 -2

(тензорного произведения) собственных векторов операторов ох и 8 г:

Щ [р] = Ух [Рх ] ® Уу [Ру ] , что в координатной записи соответствует формуле

Жи[Рх,Ру](кх,ку) = Ух[Рх](кх[Ру](ку) , (37)

Аналогично при п = 3 взаимно ортогональные собственные векторы оператора 1)3 образуют базис в указанном выше 8 • Мх• Му • М* -мерном линейном

пространстве К^ периодических 3-мерных таблиц и могут быть представлены в

л2 л2

виде тензорного произведения собственных векторов операторов Ъх, 5У и 8г:

Щ [р] = Ух [Рх ] ® Уу [Ру ] ® У* [р* ], что в координатной записи соответствует формуле

Собственные векторы Ух[рх] е К2Л ": , е К2Л-Г и [рг] е Ш2^

л2 л2 л2

операторов дх , 8У и б2 соответственно определяются, согласно (31), следующими формулами для их координат:

^а[Ра](ка) = ^(лРака/Ма) При р« = м« ,

^а[Ра](ка) = ^(ЛРака1Ма) пГи Ра = Мх + 1,...,2Ма - 1

к« = 0,1,...,2М«-1, а = х,у,* .

(39)

Рассмотрим решение уравнений (32) и (33) в периодическом случае. Наиболее важным шагом в ходе численного решения является нахождение координат -коэффициентов Фурье (обозначим их через F [p]) массива значений на сетке

правой части уравнения Fh = )} в базисе собственных векторов [W!t [р]\.

которые определяются формулами (37)-(41), стоящего в левой части уравнений

(32) и (33) оператора Dn , которые определяются через скалярное произведение

в пространстве по формуле, аналогичной (25):

Пр] = {ГнШр]\/(Щ[р1Щ[р] \ ■ (40)

Отметим, что скалярное произведение (Glh,G2/г в пространстве сеточных массивов с элементами Gsh (к) определяется формулой

G,G2h )h =^Glh(к)-G2h(к), к

которая в зависимости от размерности задачи n принимает вид

2 Nx-1 2 Ny-1

(G1h. G2h )h = ^ ^ Glh (kx, ky ) • G2h (kx, ky ) . n = 2 .

kx =0 ky =0

2 Nx-12 Ny -1 2 Nz-1

(G1h.G2h )h = X 2 2 Glh (kx'ky'ky°2h (kx'ky'kz)

n = 3.

kx =0 ky = 0 к = 0

(41)

В ходе численного решения

Тогда правые части уравнений (32) и (33) с компонентами

Фь(к) = Рь(к) + -\Ъх+Ъууь(к), п = 2;

Фк{к) = Рк{к) + -\Ъх+Ъу+Ъгу11{к), п = 3

в базисе определяемых формулами (37)-(39) собственных векторов {Щ [р]} оператора 11 п в левой части уравнений (32) и (33) будут иметь координаты

Ф[ p ] = [1 +1 (Ь x [Px ] + Ь У [ Py ])) F [ p ], n = 2, Ф[р ] = [1 + 71 (bx [Px ] + by [Py ] + 4 [ pz ])] F [p], n = 3.

(42)

Условие разрешимости уравнений (32) и (33) полностью аналогично условию разрешимости исходной задачи (13) и состоит в ортогональности правой части

первой собственной функции = 1 в пространстве И^ , то есть имеет

>

вид (/^[0]) =0. В координатной форме в зависимости от размерности

задачи п это условие имеет вид 2Щ-\ 2Му-1

2 0 при

=0 Ау =0

2 Nv -12 Ny -1 2 Nz -1

2 2 = Ъ при п = 3

=

=0 =0 к„ =0

Аналогично формуле (26) для непрерывной исходной задачи (13) координаты и [ р] (коэффициенты Фурье) решения уравнений (32) и (33) в базисе (37)-(39)

собственных векторов [р]} оператора )п в левой части уравнений (32) и

(33) при ||р|| > 1 определяются формулой

и[р] = 12Ф[р]Д[р]. (43)

Подстановка в нее формул (35) и (42) дает выражение для и [ р] в зависимости от размерности задачи п . Тогда, если для удобства считать, что и[0] = 0 , то сеточные компоненты решений уравнений (32) и (33) аналогично (27) определяется формулами

2Их-\ 2Му~1 2Ыг-1 Рх =0 Ру =0 Рг =0

2^-1 ^у-1

и.

и и

(44)

Рх =0 Ру =0

2Nx-\ 2Ny~l 2Nz-l Рх =0 Ру =0 Pz =0

■vx\px](kx}vy\py](ky}v2\p2](k2), n = 3.

Для получения сеточного массива решения Uh (к) на исходной сетке используем

обратную интерполяцию, которая дается формулами (17) и (18). Формулы (35)-(44) позволяют построить алгоритм решения уравнений (32) и (33) в периодическом случае, которые являются дискретным аналогом задачи (13).

5. Результаты тестирования

Мы создали набор программ на языке FORTRAN с использованием OpenMP, которые реализует описанный выше алгоритм численного решения краевых задач Дирихле, Неймана и периодической задачи в 2-х и 3-мерных прямоугольных областях, размеры которых удовлетворяют условию (2).

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

большим числом узлов порядка 106 -109 имела место загрузка всех ядер на 100% практически на все время работы программы, причем 3-мерная задача на

сетке с числом узлов 109 решается на 4-ядерном (с 8-ю нитями) персональном компьютере i7 примерно за 1 минуту.

Точность программ проверялась на задачах с известным аналитическим

решением по следующей схеме. В исходной области Пп,1 задаются по явным формулам потенциальное электрическое поле Е (х) = -У(( х) и его потенциал ((х) , а также соответствующая им правая часть /(х ) = Д((х) . В области Пп,1 задавалась исходная сетка, не удовлетворяющая условию (1), и массив правой части /(х) на этой сетке. Далее задача решалась численно с 4-м порядком точности, в результате чего вычисляется массив значений потенциала на исходной сетке ((к) . По нему при помощи формул разностного дифференцирования 4-го порядка точности в узлах сетки вычислялся массив значений электрического поля Ек (к ) = -Vи ((к) , который сравнивался с

точными значениями электрического поля в узлах сетки Е (г (к)) .

Относительная погрешность в вычислении компонент электрического поля определялась формулой

^(тАах(Еа(г(к))"ЕаиИк)))/тах(Е«(г(к))) .

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

((г (к)) при помощи формул 4-го порядка точности, которая определяется

формулой

И ^па^(тах(Е«(г(кг(к))))/тах(Е«(г(к))) .

Многочисленные тесты показали, что если сетка достаточна для отображения изменения потенциала и его лапласиана, то эти погрешности имеют одинаковый

порядок, и погрешность превосходит погрешность 6и не более чем в 2 раза,

то есть верны оценки

Ек < 2, ^ - 1.

¿и

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

= тах

а=х,

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

Литература

1. Владимиров В.С. Уравнения математической физики. 5-е изд-е. Москва: Наука,

1988.

2. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

3. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва: Лаборатория базовых знаний, 2001, 632 с.

4. Деммель Дж. Вычислительная линейная алгебра. Теория и приложения.- М.: Мир, 2001.

Сведения об авторах Мингалев Олег Викторович

к.ф.-м.н., зав.сектором, Полярный геофизический институт, Апатиты E-mail: mingalev_o@pgia.ru

Мельник Михаил Николаевич

м. н. с., Полярный геофизический институт, Апатиты E-mail: melnik@pgia.ru