АЛГОРИТМЫ И ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ
УДК 514.75(08)
С. А. Ишанов, С. В. Клевцур, А. И. Кожурова, К. С. Латышев, В. Н. Худенко
ЦИКЛИЧЕСКИЙ ВАРИАНТ «а-р» ИТЕРАЦИОННОГО МЕТОДА. ОЦЕНКИ СКОРОСТИ СХОДИМОСТИ АЛГОРИТМА
ного уравнения диффузии ионов со смешанными производными и первыми производными дивергентного вида. Проведены тестовые расчеты на модельной задаче с известным аналитическим решением. Показана работоспособность алгоритма и дана оценка скорости его сходимости.
Differences methods of the approximate solution of the two-dimensional equation of ions diffusion with the mixed derivatives and the first derivatives of a divergent look are considered. Test calculations on a modeling task with the known analytical decision are carried out. Operability of algorithm is shown and the estimation of speed of its convergence is given.
Ключевые слова: разностный метод, приближеное решение, уравнение с частными производными, «а — Р» алгоритм, сходимость.
Key words: differences method, approximate solution, equation with private derivatives, «а — Р» algorithm, convergence.
При постановке задачи глобального моделирования ионосферы в сферической географической системе координат в качестве граничных условий по одной или нескольким пространственным переменным могут быть использованы условия периодичности решений и коэффициентов исходной системы уравнений. Возникает необходимость в разработке численных алгоритмов, учитывающих такие особенности задачи.
Рассмотрим уравнения диффузии ионов О+ и Н+, записанные в дивергентной форме в полярной системе координат:
где N1 — концентрация ионов и скорость движения ионов О+ и Н+ (і = 1, 2); Qi, и Li — скорость образования и коэффициент, характеризующий рекомбинационные процессы; t — время; г — координата вдоль радиуса-вектора; X — полярный угол; Prr,..., PXr — коэффициенты дифференциального оператора, относящиеся к иону сорта і. Выражения для этих коэффициентов приведены в работе [1].
Используя разностные схемы из [2 — 3], уравнение (1) может быть представлено в виде системы девятиточечных разностных уравнений
68
Рассмотрены разностные методы приближенного решения двумер
(1)
Вестник Балтийского федерального университета им. И. Канта. 2012. Вып. 10. С. 68 —72.
^Уі,к Аі,куі-1,к-1 + Ві,куі,к-1 + LifkУi+1fk-1 + Кі,куі-1,к (2)
— Сі,куі,к + Еі,куі+1,к + Ц,куі-1,к+1 + Vilkyilk+1 =—Fifk,
2^гЧЛ/, -1Д = 0,+1,+2,....
Отметим, что коэффициенты и правая часть разностных уравнений (2) периодичны по индексу к с периодом N.
Одним из наиболее эффективных способов решения системы разностных уравнений (2) является итерационный «а— р» алгоритм [3], который представляет собой двумерный вариант метода прогонки по двум пространственным координатам. Данный метод сохраняет устойчивость при резких и зачастую трудно предсказуемых изменениях коэффициентов разностных уравнений (2), не требует получения какой-либо априорной информации о границах спектра разностного оператора, что необходимо для эффективного применения итерационных других методов [4], и имеет высокую скорость сходимости.
В то же время «а — р» ■итерационный метод не может быть использован непосредственно для решения двумерных разностных задач с периодическими краевыми условиями. В связи с этим был разработан модифицированный вариант этого алгоритма для решения периодических краевых задач.
Дополним систему разностных уравнений (2) периодическими граничными условиями:
У,,и = ак, УК д =Ьк, і ак = ак+щ, Ък = Ък+Щ. (3)
В этом случае решение системы (2), если оно существует, тоже будет периодическим с периодом N, то есть
У и = У.д+М,' 1 ^ М'. (4)
Поэтому достаточно найти решение уікг например при к = 1, 2, ..., М' -1, учитывая, что у)Л = у, ... ,1 $; < Л/(.
Так же как и в работе [4], в которой рассмотрен одномерный случай, решение будем искать в виде линейной комбинации сеточных функций Хі,к и 2ід:
У і,к = х,-,* +1/, д=,д, 1 ї' і < Л/„ 1 < к .< Ык. (5)
Функцию хі,к определим как решение неоднородной системы разностных уравнений
Лха = -Р1:к, 2 ^ -1, 2 ^ ЛГ* -1 (6)
с однородными условиями по периодичности
х,,1 = *і,щ = 0, х,д = йк, хК: д = Ък, 2 ^ гЧ Ц -1,1 ^ к ^ Ык . -1. (7)
Функцию zi к определим как решение однородной системы разностных уравнений
Лга = 0, 2 ^ гЧ Ц -1, 2 «с к^ -1 (8)
с неоднородными условиями по периодичности
~і,і = ~іж = 1' -'ід = ,к = 2 ^ і ^ М/ — 1,1 ^ к ^ Ык — 1. (9)
69
70
Запишем граничные условия для системы (6) в общем виде, тогда с учетом (7) получим:
ф2,к = °, Ф і,2 = О, Ф N-1,к = °, ф1 N -1 = °,
' 1 (10)
£>2,к = ак, -1,к = Ьк, ^1,2 = 0, -1 = °.
Аналогично, для системы (8) с учетом (9) находим:
ф2,к = °, фЫ, -1,к = °, ф1,2 = °, ф1 ,К, -1 = 0
(11)
^2,к = °, ^ -1,к = °, §і,2 = 1, ,Ык-1 = 1.
Пусть решение системы (6) — (7) удовлетворяет соотношениям: хі ,к = а1+1, кХі+1,к + р+1,к , хі ,к = У і-1, кХі-1,к + йі-1,к , хі,к = аі,к+1Хі,к+1 +рі,к+:^ хі,к = Уі,к-1 Хі ,к-1 + йі ,к-1.
Для отыскания неизвестных прогоночных коэффициентов а, у, р, й, а, ~ , р , й применим «а—р» итерационный алгоритм, сводящийся к решению восьми нелинейных алгебраических уравнений для каждого узла разностной сетки [5]. Определив пару значений прогоночных коэффициентов, легко найти решение задачи (6) — (7). Аналогично определяется решение задачи (8) — (9).
Отметим, что прогоночные коэффициенты а, у, а, у для задачи (8) — (9) рассчитываются по тем же формулам, что и в задаче (6) — (7), и при тех же граничных условиях (1°) — (11). Таким образом, общее число «а — р» уравнений для прогоночных коэффициентов в двух данных задачах сокращается до двенадцати.
Легко показать, что решение у,,к, определенное по формуле (5), удовлетворяет исходной задаче (2) — (4) во всех узлах разностной сетки, кроме узлов С номером к = I. Определим решение у і Д (1 '' і < Ы,), для
этого подставим решение (4) в уравнение (2) при к = 1 и приведем полученное выражение к трехточечному виду относительно уг 1:
А,1уг-1,1 - сг ду1,1 + ду1+1,1 = -~ ,1, (12)
где
А,1 = К,1 + А,121 -1,Nk-1 + Д,121 -1,2 , Ci,1 = С1,1 - В1,121 N -1 - ^1,121,2 ,
~1,1 = Е ,1 + Ц ,12г+1,Nk -1 + ^ ,12г+1,2'
^ ,1 = ^,1 + А ,А-1,Nk-1 + ,1Хг ,Nk-1 + ,1Хг+l,Nk-1 +
+ Ц ,1Х1 -1,2 + ,1Хг-1,2 + ,1Хг+1,2"
Видно, что коэффициенты и правая часть трехточечного уравнения
(12) могут быть вычислены по ранее найденным значениям сеточных функций хг к и к. Так как уравнение (12) одномерно, то для отыскания его решения можно воспользоваться формулами обыкновенной прогонки [4] с граничными условиями
уг,1 = Й1, yN, ,1 = Ь1. (13)
= с=а*
Таким образом, циклический вариант «а — р» итерационного алгоритма состоит из трех этапов:
1) методом «а— р» итераций находим значения сеточных функций хг к и ,к во всех внутренних углах сеточной области;
2) находим значения функции уц (1 ?: і Л/,), решая одномерной
прогонкой уравнение (12) с граничными условиями (13);
3) по найденным значениям сеточных функций х, к и 2,д находим значения уі к во всех внутренних узлах сеточной области по (5).
Для проверки работоспособности и оценки скорости сходимости модифицированного алгоритма рассмотрим тестовую задачу с известным аналитическим решением для эллиптического уравнения, записанного в полярной системе координат. В декартовой системе координат область определения решения — кольцо с радиусами Я° и Яі. В полярной системе координат области соответствует область Є = {Я° < г < Я1,11 < X < 12,11 > °, 12 -11 = 2я|. Найдем решение уравнения
1 ^ диІ 1_ д^и г дг ^ дг ) г2 дХ2
Г1ГІ+-= Кг, X). (14)
Решение является периодическим по X с периодам 2л и удовлетворяет на границе области О краевым условиям первого рода:
и(г, X) = й(Х), г = Яд, X е [11,12), и(г, X) = Ь(Х), г = Я1, X е [11,12).
Правая часть / (г, X) подбиралась так, чтобы функция
и(г, X) = (г - Я0 )(г - Я1 )8т X+С,
где С — некоторая константа, удовлетворяла уравнению (14), при этом а = Ь = с. Постоянная С выбиралась в пределах 0 ь 105. Разностная аппроксимация уравнения (14) проводилась по схемам из [2 — 3].
Вычислительные эксперименты показали хорошее совпадение численного и точного решений. Погрешность расчетов не превосходила 10%. В таблице дана зависимость числа «р» итераций, необходимых для достижения точности е = 10 3, от числа узлов разностной сетки (число «а» итераций во всех расчетах не превосходила 4). Здесь Ni — количество узлов по г, N — по X.
Зависимость числа «р» итераций от числа узлов разностной сетки
71
Число узлов 1° 1° 2° 2° 3° 3° 3° 5°
N 19 73 19 37 19 37 73 73
Число «р» итераций по х,, к 15 29 22 39 23 53 1°9 67
Суммарное число «р» итераций 26 55 37 7° 4° 91 192 344
Суммарное число «р» итераций получается сложением числа «р» итераций, затраченных на вычисление прогоночных коэффициентов Х, к и 2, к . Все столбцы таблицы 1, кроме последнего, где с = 1°5, соот-
72
ветствуют с = °. Анализ полученных результатов показывает, что скорость сходимости циклического «а - р» итерационного метода остается достаточно высокой. При этом не требуется ни самосопряженности разностного оператора, ни априорной информации о его спектре.
Скорость сходимости асимптотически уменьшается с ростом числа узлов сетки пропорционально I / • Ык, как и в исходном алгоритме.
Циклический алгоритм решения двумерных периодических задач ионосферного моделирования не является жестко связанным с той или иной конкретной методикой расчета уравнений и может быть использован практически без изменений в целом ряде программных комплексов, в которых для расчета используются четырехугольные сетки.
Работа выполнена при финансовой поддержке РФФИ по проектам № 11-01-00098а и № 11-01-00558а.
Список литературы
1. Фаткуллин М. Н., Клевцур С. В., Латышев К. С. Оператор переноса в уравнении непрерывности для ионов в трехмерно-неоднородной области F (средние и высокие широты) / / Геомагнетизм и аэрономия. 19S4. Т. 24, № б. С. 90б — 910.
2. Ишанов С. А., Клевцур С. В, Латышев К. С. Алгоритм «а—p» итераций в задачах моделирования ионосферной плазмы // Математическое моделирование. 2009. Т. 21, № 1. С. 33—45.
3. Ишанов С. А., Клевцур С. В. Математическое моделирование ионосферы с учетом ее трехмерной неоднородности // Вестник Российского государственного университета им. И. Канта. 2010. № 4. С. 152 — 15S.
4. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М., 197S.
5. Четверушкин Б. Н. Математическое моделирование задач динамики излучающего газа. М., 19S5.
Об авторах
Сергей Александрович Ишанов — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта.
E-mail: sergey.ishanov@ya.ru.
Сергей Владимирович Клевцур — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта.
Алла Ивановна Кожурова — ст. преп., Балтийский федеральный университет им. И. Канта.
Константин Сергеевич Латышев — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта.
Худенко Владимир Николаевич — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта.
Authors
Dr Sergey Ishanov — professor, I. Kant Baltic Federal University.
E-mail: sergey.ishanov@ya.ru.
Dr Sergey Klevtsur — assistant professor, I. Kant Baltic Federal University.
Alla Kozhurova — high instructor, I. Kant Baltic Federal University.
Dr Konstantin Latyshev — professor, I. Kant Baltic Federal University.
Dr Vladimir Khudenko — assistant professor, I. Kant Baltic Federal University.