CC BY
61
Вычислительные технологии
Том 16, № 1, 2011
Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнения Пуассона*
В. И. Исаев1, В. П. Шапеев2, С. В. Идимешев1 1 Новосибирский государственный университет, Россия 2Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected], [email protected]
Предложен способ построения вариантов метода коллокаций и наименьших квадратов высокого порядка точности для численного решения уравнения Пуассона, на основе которого созданы варианты метода до восьмого порядка точности включительно. В численных экспериментах показано, что в случае достаточно гладкого решения эти варианты позволяют получать приближенное решение, сходящееся к точному с высоким порядком на последовательности сеток при h ^ 0, где h — максимальный линейный размер ячеек сетки.
Ключевые слова: численные методы, метод коллокаций и наименьших квадратов, повышенный порядок точности, уравнение Пуассона.
Введение
Проблема создания методов повышенного порядка точности численного решения эллиптических уравнений успешно решалась многими исследователями начиная с 1930-х годов [1-6]. В данной работе впервые предлагаются варианты метода коллокаций и наименьших квадратов (КНК) повышенного порядка точности для численного решения уравнения Пуассона. Это уравнение как достаточно простое было выбрано для отработки подхода к созданию вариантов метода КНК повышенного порядка в более сложных случаях.
Метод коллокаций успешно применялся для численного решения краевых задач для обыкновенных дифференциальных уравнений (ОДУ) [7-11] и уравнений с частными производными [12-14]. Суть метода заключается в следующем. Приближенное решение ищется в конечномерном линейном пространстве функций. Неизвестные коэффициенты его разложения по базису пространства определяются из уравнений коллокаций и краевых условий. Уравнения коллокаций представляют собой требования того, чтобы приближенное решение удовлетворяло уравнениям исходной дифференциальной задачи в конечном множестве точек области (в точках коллокаций), в которой ставится данная задача. Краевые условия получаются из соответствующих условий задачи, записанных в нескольких точках на границе области. В методе коллокаций рассматривается ровно столько уравнений, сколько имеется неизвестных. В методе КНК число уравнений
*Работа выполнена при финансовой поддержке РФФИ (гранты № 08-08-00249-а и 10-01-00575-а), Интеграционных проектов СО РАН № 11.5 и 26 и Комплексного интеграционного проекта СО РАН № 140.
превышает количество неизвестных, т, е, система, из которой вычисляются неизвестные коэффициенты, является переопределенной. Для ее решения используется метод наименьших квадратов (\ 111К).
В [13,15] предложены и реализованы варианты методов коллокаций и КНК для эллиптических уравнений второго порядка, В них приближенное решение ищется в пространстве кусочно-квадратичных функций на сетках с прямоугольными или треугольными ячейками, В настоящей работе на основе результатов, полученных в [12,13,15-18], предложен способ построения вариантов метода КНК высокого порядка точности для уравнения Пуассона, на основе которого созданы варианты метода до восьмого порядка точности включительно. Численные эксперименты показали, что при достаточно гладком решении применение этих вариантов позволяет получать приближенное решение, сходящееся к точному с высоким порядком на последовательности сеток при к ^ 0, где к — максимальный линейный размер ячеек сетки,
1. Постановка задачи и описание метода
В области О = [0,1] х [0,1] рассмотрим задачу Дирихле для уравнения Пуассона
д 2 и д 2 и
Щ + Щ = ^ (хг,х2)еП, (1)
и = д, (х1,х2) Е дО*, (2)
где f (х1 ,х2) и д(х1,х2) — известные функции, а и(х1,х2) — искомое решение. На примере решения этой задачи будем рассматривать алгоритм метода КНК, В случае других краевых задач для уравнения Пуассона реализация метода осуществляется аналогичным образом.
Покроем область О сеткой, состоящей из N прямоугольных ячеек О1, ..., Ом. Пусть 2кн, 2к2г — размеры яч ейки О, то направлениям осей х1 и х2 соответствен но, I = 1, ..., N. Метод КНК используется в данной работе в сочетании с методом декомпозиции области [19-21], В последнем решение задачи (1), (2) сводится к построению сходящегося итерационного процесса, на каждом шаге которого требуется последовательно решать вспомогательные локальные задачи в подобластях О1, ..., Ом [19-21], Отметим, что в методе декомпозиции предлагается естественный способ распараллеливания вычислительного алгоритма, поскольку в ходе итерационного процесса данные обрабатываются раздельно по подобластям [21, 22], Локальная задача в ячейке О, имеет вид
д2ик д2ик
+ = (жьж2)еа, (3)
дик ди
Пк + кг—=й + кг —, {Хг,Х2) Е дПг \ Ш, (4)
дп дп
ик = д, (х1 ,х2) Е дО, П дО, (5)
где ик — приближенное решение, полу чаемое на итерации, и — приближение для и, имеющееся в соседних ячейках при построении решения на к-й итерации в подобласти О,, п = (п1, п2) — вектор единичной внешней нормали к границе дО,, к, = (к1,к2,)1/'2,
дпк дпк дпк
т— = п 1---Ь п2т;—•
ОН дх1 дх2
Условие (4) записывается в данной работе по аналогии с предложенным в [15],
Для удобства реализации численного метода введем в ячейке Qi локальную систему координат
>= Х1 х1сг х2 Х2Ы / п\
- = —КГ' 6 = (6)
где (х1Сг,Х2сг) _ центр ячейки, После замены переменных уравнения (3)-(5) примут вид
д2 ьк д2 ьк ~ -
^ + = (6,6) е а, (7)
дьк „ дь
где
+ ^ = v + (СьЫеШЛШ, (8)
vk = g, (6,6) е дйг n дП, (9)
vk (6,6) = uk (xi (6,6)x (6,6)) , v(Cl,C2) = U (xi (6,6) x (6 ,6)) , fiCbb) = f{x1(£l,b),x2 ^ъЫ) , ЖъЫ = g(x1 ^ъб^^б^) ,
dvk dvk dvk ni^r + n2-
hu hu X h2i h2i~
hi hi hi hi
dn d^l д£2
Здесь Qi, Q — образы множеств Qi и Q при аффинном отображении, задаваемом формулами перехода к локальной системе координат (6),
Qi
Для численного решения локальных задач (7) (9) применяется метод KHK, В предлагаемых здесь вариантах метода в каждой ячейке сетки Ql; ..., QN приближенное решение ищется в пространстве Pm(Qi) полиномов степени не выше m, определенных Qi
СГ СТ, 0 < al < m, 0 < а2 < m - al. (10)
Всего в каждой ячейке сетки имеем d = (m+l)(m+2)/2 базисных элементов. Обозначим их pl, ..., pd. Приближенное решение в ячейке Qi представляется в виде линейной комбинации
d
v = ^ Cik Pk, k=l
где Cil, ..., Cid — искомые коэффициенты разложения по базису, которые находятся из переопределенной системы линейных алгебраических уравнений (СЛАУ), состоящей из уравнений коллокаций и условий на границе ячейки. Уравнения коллокаций здесь представляют собой требования того, чтобы приближенное решение удовлетворяло (7)
Qi
дQi — требования выполнения уравнений (8), (9) в конечном числе узлов. Среди этих
требований можно выделить те, которые записываются па границах между соседними ячейками (8). Будем называть их условиями согласования (условиями сшивки |21|). В случае д Пг П дП = 0 в нескольких точках пере сечения дПг П д П записываются краевые условия из исходной постановки задачи (2),
В данной работе точки записи уравнений коллокаций, условий согласования и краевых условий распределялись равномерно внутри каждой ячейки и на ее границе (рис, 1), Их количество выбиралось таким образом, чтобы в СЛАУ дня определения коэффициентов разложения но базису линейно независимых уравнений коллокаций было не менее т (т — 1)/2, а независимых условий на границе — не менее 2 т +1. Здесь т (т — 1)/2 и 2 т +1 — соответственно размерности образа и ядра оператора Лапласа, отображающего Рт(Пг) на Рт-2(Пг). Общее число независимых уравнений в СЛАУ в предложенных вариантах метода КНК превышает количество неизвестных. Здесь были реализованы варианты с различными соотношениями числа уравнений и количества неизвестных. Дня всех вариантов, используемых ниже в раздело 2,1, числа точек коллокаций, точек записи условий согласования и краевых условий, а также их расположение совпадают с указанными на рис, 1,
Дня решения переопределенной СЛАУ вместо МНК используется ортогональный метод линейной алгебры |23|, который дает то же решение, что и МНК при отсутствии ошибок округления. Преимущество ортогонального метода заключается в следующем. Пусть рассматривается переопределенная система
где А — прямоугольная вещественная матрица, х, Ь векторы неизвестных и правых частей соответственно, В методе наименьших квадратов ее решение находится из условия минимума функционала невязки
Здесь и далее (а, Ь) — скалярное произведение векторов а и Ь, В известном варианте реализации МНК |24| фактически решается СЛАУ
А х = Ь,
(Н)
Ф(х) = (Ах — Ь, Ах — Ь).
АТА х = Ат Ь,
(12)
Рис. 1. Схема расположения точек коллокаций (о), точек записи условий согласования (□) и краевых условий (•) в ячейке сетки
здесь AT — транспонированная матрица A. Заметим, что
* (ATA) = (v (A))2,
где V(A) = amax/omvn — число обусловленности матрицы A, amax, amin — ее максимальное и минимальное сингулярные числа. Таким образом, обусловленность задачи (12) может быть существенно хуже обусловленности системы (11). В ортогональном методе СЛАУ (11) приводится к виду
R x = Q b, (13)
где R = QA — верхнетреугольная матрица, т. е. ее элементы, расположенные ниже глав-
Q
(11) и (13), полученные методом наименьших квадратов, совпадают при отсутствии ошибок округления. Для (13) это решение можно легко найти из системы с верхнетреугольной квадратной матрицей, состоящей из первых n уравнений СЛАУ (13), где n — размерность вектора x. Заметим, что v(R) = v(A), так как R = QA. Таким образом, ортогональный метод в отличие от МНК не ухудшает обусловленность системы уравнений в процессе ее решения.
Для уменьшения времени, необходимого для проведения расчета, в данной работе использовался метод ускорения сходимости итераций, в котором в ходе итерационного процесса к текущему приближению добавлялась поправка из подпространства Крылова [25-28].
2. Численные эксперименты
2.1. Сходимость приближенного решения на последовательности сеток
Для исследования порядка сходимости приближенного решения uh(xi,x2) к точному проведены расчеты на последовательности сеток при h ^ 0 где h — максимальный линейный размер ячеек сетки. Ниже приведены результаты численных экспериментов с решением
u(xi,x2) = exp(10 xi) + exp(10 x2) + exp(5 xi + 5 x2), (14)
которое удовлетворяет уравнению Пуассона (1) с правой частью
f (xi ,x2) = 100(exp(10 xi) + exp(10 x2)) + 50exp(5 xi + 5 x2).
Решение (14) в области Q имеет весьма большие градиенты. Здесь оно используется для демонстрации возможностей вариантов метода KHK, предложенных в данной работе. Расчеты проводились на равномерных регулярных сетках с ячейками
Qij = [(г - 1) 2hi, г 2hi] х [(j - 1) 2h2, j 2h2], г =1,...,Ni, j = 1, ..., N2,
где Nl и N2 — размеры сетки по направлениям осей xi и x2 соответственно, 2hi = 1/Ni, 2h2 = 1/N2. В таблице приведены значения погрешности
eni,N2 = \\u - uh\\х = sup |u(xi,x2) - Uh(xi, x2)|
(х1,Ж2)€П
приближенных решений, построенных вариантами метода KHK при m = 2, ..., 8. Для наглядности приведены отношения погрешностей, полученных на двух последовательных сетках размера Ni х N2 и Ni/2 х N2/2. Для m = 2 погрешность ENltN2 на сетке
Погрешности приближенных решений, полученных методом КНК
т N1 X N2 Ей !,м2 ЕN-1/2, N2/2 ЕМг ,М2
2 80 х 80 160 х 160 320 х 320 20.7 5.18 1.29 4.0 4.0
3 80 х 80 160 х 160 320 х 320 2.95 7.36 • КГ1 1.84 • Ю-1 4.0 4.0
4 20 х 20 40 х 40 80 х 80 5.20 • Ю-1 3.26 • 10"2 2.04 • 10"3 15.9 16.0
5 20 х 20 40 х 40 80 х 80 1.19 • Ю-1 7.39 • 10"3 4.61 • 10"4 16.0 16.0
6 10 х 10 20 х 20 40 х 40 5.52 • 10"2 8.98 • 10"4 1.41 • 10"5 61.4 63.7
7 10 х 10 20 х 20 40 х 40 9.81 • 10"3 1.53 • 10"4 2.37- 10"6 64.2 64.4
8 10 х 10 20 х 20 40 х 40 2.43 • 10"4 9.89 • 10"7 4.41 • 10"9 245 224
80 х 80 равна 20.7. При этом относительная погрешность здесь составляет менее 0.04%. Видно, что для рассматриваемого гладкого решения с ростом т существенно повышается точность расчета. Приведенные в таблице результаты свидетельствуют о том, что порядок сходимости приближенного решения к точному при к ^ 0 здесь не хуже т-го при четных т и (т — 1)-го при нечетных т.
2.2. Влияние погрешностей округления на приближенное решение
Известно, что погрешность приближенного решения на ЭВМ математической задачи обусловливается несколькими факторами [29]: погрешностью метода, вычислительной погрешностью, вызванной ошибками округления при арифметических действиях, неустранимой погрешностью, являющейся следствием неточности входных данных. С ростом числа ячеек сетки уменьшается погрешность сеточного метода. Одновременно с этим увеличивается число арифметических действий, выполняемых при построении решения, и происходит накопление погрешностей округления. При фиксированной разрядности представления вещественных чисел на ЭВМ наилучшая точность приближенного решения, которая может быть достигнута при помощи выбранного численного метода, зависит от устойчивости решения к возмущениям, возникающим из-за неточности входных данных и ошибок округления.
Рассмотрим модельную краевую задачу для ОДУ второго порядка
— =-104 8т(100ж), же (0,1), и{0) = 0, м(1) = вт(100).
(15)
На ее примере проведем сравнение приближенных решений, получаемых методами кол-локаций и КНК, Варианты этих методов дня задачи (15) реализовывалиеь здесь так же, как и в разделе 1, В них приближенное решение локальных задач, аналогичных (7)-(9), в каждой ячейке сетки ищется в пространстве Pm([— 1,1]). В одномерном случае ячейка является отрезком. Ее граница состоит из двух точек, В зависимости от расположения ячейки в этих точках записываются условия согласования или краевые условия задачи (15),
Рассмотрим два варианта метода коллокаций. В первом в качестве точек коллокаций в каждой ячейке сетки использованы корпи полиномов Чебышева первого рода
( 21 — 1 \ 6 = 008 Г 2(^1))' * =
где £ — локальная переменная в ячейке, которая определяется аналогично (£ь£2) в (6)-Во втором варианте точки коллокаций распределены равномерно в каждой ячейке |8|
t 21 ~m 1 1 1 6 =-7Г, / = 1, .. •, те — 1.
m — 2
Дня сравнения рассмотрим вариант метода КНК, в котором точки коллокаций также распределены равномерно:
а = --1, / = 0, ...,2т.
m
При этом число уравнений, записываемых в каждой ячейке, примерно в 2 раза превышает число неизвестных.
На рис. 2 приведены графики зависимости погрешности
En = \\u — uh= sup \u(x) — uh(x)\
0<x< 1
приближенных решений Uh, полученных указанными методами при m = 20, от числа ячеек сетки N. Для наглядности по обеим осям используется логарифмический масштаб. Начиная с некоторого размера сетки погрешность (см. рис. 2) с ростом числа ячеек перестает убывать. Это происходит, когда возмущения приближенного решения,
Рис. 2. Зависимость от числа ячеек сетки N погрешности приближенных решений, полученных методом КНК (1) и вариантами метода коллокаций, в которых точки коллокаций расположены равномерно (второй вариант, 2) и в корнях полиномов Чебышева (первый вариант, 3) в каждой ячейке сетки
возникающие из-за неточности входных данных и ошибок округления, становятся соизмеримыми с погрешностью метода, В данной работе для представления чисел на ЭВМ использовался тип данных double (С++), в котором мантисса состоит из 15 разрядов. Из графиков на рис, 2 следует, что наиболее точные приближенные решения, которые удается получить рассматриваемым здесь вариантом метода КИК и первым вариантом метода коллокаций, имеют примерно одинаковую погрешность EN & 10_13, Это, по-видимому, свидетельствует о том, что погрешности округления в методе КИК и в первом варианте метода коллокаций влияют на точность приближенного решения примерно в равной степени. Для второго варианта метода коллокаций минимальное значение погрешности EN & 2 • 10_11, Отсюда можно предположить, что во втором варианте метода коллокаций погрешности округления сильнее влияют на точность решения, чем в первом варианте и в методе КИК, Отметим, что аналогичная ситуация наблюдается в задаче аппроксимации функции полиномами [30,31], Известно, что при использовании корней полиномов Чебышева в качестве узлов интерполяции устойчивость полинома Лагранжа к погрешностям входных данных и возмущениям, возникающим из-за ошибок округления, существенно лучше по сравнению со случаем использования равномерно распределенных (равноотстоящих) узлов [30, 31], Применение метода наименьших квадратов для построения аппрокеиманта также, как правило, позволяет добиться его хорошей устойчивости. Это, несомненно, связано с тем, что в МНК аппрокеимант находится из условия минимума функционала невязки.
Список литературы
[1] Микеладзе Ш.Е. О численном интегрировании уравнений эллиптического и параболического типов // Изв. АН СССР. Математика. 1941. Т. 5, № 1. С. 57-74.
[2] Валиуллин А.Н. Схемы повышенной точности для задач математической физики. Новосибирск: НГУ, 1973.
[3] Валиуллин А.Н., Ганжа В.Г., Ильин В.П. и др. Задача автоматического построения и исследования на ЭВМ разностных схем в аналитическом виде // Докл. АН СССР. 1984. Т. 275, № 3. С. 528-532.
[4] Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск: Изд-во Ин-та математики СО РАН, 2000.
[5] Шапеев A.B., Шапеев В.П. Разностные схемы повышенной точности для решения эллиптических уравнений в области с криволинейной границей // Журн. вычисл. математики и мат. физики. 2000. Т. 40, № 2. С. 223-232.
[6] Шапеев В.П., Шапеев A.B. Решение эллиптических задач с особенностями по схемам высокого порядка аппроксимации // Вычисл. технологии. 2006. Т. 11. Спец. выпуск, посвященный 85-летию со дня рождения акад. H.H. Яненко. Ч. 2. С. 84-91.
[7] Russell R.D., Shampine L.F. A collocation method for boundary value problems // Numer. Math. 1972. Vol. 19. P. 1-28.
[8] de Boor C., Swartz B. Collocation at Gaussian points // SI AM J. Numer. Anal. 1973. Vol. 10, No. 4. P. 582-606.
[9] Ascher U., Christiansen J., Rüssel R.D. A collocation solver for mixed order systems of boundary value problems // Math. Comput. 1979. Vol. 33. P. 659-679.
[10] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980.
[11] Schild К.Н. Gaussian collocation via defect correction // Numer. Math. 1990. Vol. 58. P. 369-386.
[12] Плясунова А.В., Слепцов А.Г. Коллокационно-сеточный метод решения нелинейных параболических уравнений на подвижных сетках // Моделирование в механике. 1987. Т. 1(18), № 4. С. 116-137.
[13] Слепцов А.Г. Коллокационпо-сеточное решение эллиптических краевых задач // Там же. 1991. Т. 5(22), № 2. С. 101-126.
[14] Leyk Z. A Co-collocation-like method for elliptic equations on rectangular regions // Australian Math. Soc. B. Appl. Math. 1997. Vol. 38. P. 368-387.
[15] Слепцов А.Г., Шокин Ю.И. Адаптивный проекционно-сеточный метод для эллиптических задач // Журн. вычисл. математики и мат. физики. 1997. Т. 37, № 5. С. 572-586.
[16] Беляев В.В., Шапеев В.П. Метод коллокаций и наименьших квадратов на адаптивных сетках в области с криволинейной границей // Вычисл. технологии. 2000. Т. 5, № 4. С. 13-21.
[17] Shapeev V.P., Semin L.G., Belyaev V.V. The collocation and least squares method for numerical solution of Xavicr Stokes equations // Proc. of the Steklov Inst, of Math. Suppl. 2. 2003. P. SI 15 SI37.
[18] Исаев В.И., Шапеев В.П., Еремин С.А. Исследование свойств метода коллокации и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье^Стокса // Вычисл. технологии. 2007. Т. 12, № 3. С. 53-70.
[19] годунов С.к. Уравнения математической физики. м.: Наука, 1979.
[20] Мацокин A.M., Непомнящих С.В. Метод альтернирования Шварца в подпространстве // Изв. вузов. Математика. 1985. № 10. С. 61-66.
[21] Василевский Ю.В., Ольшанский М.А. Краткий курс по многосеточным методам и методам декомпозиции области. М.: МГУ, 2007.
[22] Issaev V.I., Shapeev V.P. Development of the collocations and least squares method // Proc. of the Steklov Inst, of Math. Suppl. 1. 2008. P. S87-S106.
[23] Беклемишев Д.В. Дополнительные главы линейной алгебры. М.: Наука, 1983.
[24] Форсайт Дж., Молер К. Численное решение систем линейных алгебраических уравнений. М.: Мир, 1969.
[25] Слепцов А.Г. Об ускорении сходимости линейных итераций // Моделирование в механике. 1989. Т. 3(20), № 3. С. 132-147.
[26] Saad У. Numerical Methods for Targe Eigenvalue Problems. Manchester Univ. Press, 1991.
[27] Karamyshev V.B., Kovenya V.M., Sleptsov A.G., Cherny S.G. Variational method of accelerating linear iterations and its applications // Comput. k, Fluids. 1996. Vol. 25, No. 5. P. 467-484.
[28] Исаев В.И. Алгоритм ускорения сходимости в методе коллокаций и наименьших квадратов // Вычисл. технологии. 2008. Т. 13. Спец. выпуск 4. С. 41-46.
[29] Бахвалов И.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987.
[30] Бабенко К.И. Основы численного анализа. М.: Ижевск: НИЦ "Регулярная и хаотическая динамика", 2002.
[31] Ильин В.П. Численный анализ. Ч. 1. Новосибирск: Изд. ИВМиМГ СО РАН, 2004.
[32] Jiang Bo-Nan. On the least squares method // Comput. Methods Appl. Mech. Engrg. 1998. Vol. 152. P. 239-257.
Поступила в редакцию 5 апреля 2010 г.
