close

Вход

Забыли?

вход по аккаунту

?

Многосеточные итерационные алгоритмы в методе конечных элементов с учетом численного интегрирования.

код для вставкиСкачать
Вычислительные технологии
Том 14, № 1, 2009
Многосеточные итерационные алгоритмы
в методе конечных элементов
с учетом численного интегрирования∗
Л. В. Гилева
Институт вычислительного моделирования СО РАН, Красноярск, Россия
e-mail:gileva@icm.krasn.ru
Рассматривается сеточная задача, полученная дискретизацией эллиптического
уравнения второго порядка с помощью кусочно-линейных элементов на треугольниках с использованием численного интегрирования. Для ее решения на последовательности вложенных сеток используются полный многосеточный алгоритм
на основе W -цикла и каскадный алгоритм, являющийся наиболее простой версией многосеточных методов. Дискретные задачи на более грубых сетках строятся
таким образом, что матрицы систем уравнений на соседних сетках связаны через
операторы интерполяции и проектирования. Доказано, что для обоих алгоритмов
число арифметических операций, приходящихся на одно неизвестное, для определения приближенного решения с точностью, совпадающей по порядку с погрешностью дискретной задачи с учетом численного интегрирования, не зависит от числа
неизвестных и количества сеток.
Ключевые слова: метод конечных элементов, многосеточный итерационный алгоритм, каскадный алгоритм, квадратурная формула, оценка числа операций.
Введение
Наиболее простая версия многосеточных методов — каскадный алгоритм. Он состоит
в применении традиционных итерационных методов на последовательности вложенных сеток для достижения хорошего начального приближения на каждой сетке за счет
интерполяции приближенного решения с предыдущей, более грубой, сетки. Эффективность такого подхода изучалась в работе В.П. Ильина и В.М. Свешникова [1], однако
систематические исследования на эту тему начались только в 90-х годах. Свое название
метод получил в статье П. Дойфльхарда [2], где продемонстрирована вычислительная
эффективность алгоритма. Теоретическое обоснование дано в работах [3, 4], а затем
повторено в [5].
Последующие вычислительные эксперименты и теоретические исследования продемонстрировали экономичность каскадного алгоритма для двумерных эллиптических
краевых задач, в том числе квазилинейных [6], с криволинейной границей [7] и в областях с угловыми точками [8], где приходится сгущать сетку для компенсации потери
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований
(грант № 08-01-00621а).
c Институт вычислительных технологий Сибирского отделения Российской академии наук, 2009.
°
∗
34
35
Л. В. Гилева
гладкости. В трехмерных краевых задачах алгоритм оказался даже более эффективным
[9, 10].
В перечисленных работах сходимость каскадного алгоритма обосновывалась в предположении, что элементы матрицы и вектора правой части системы уравнений метода конечных элементов вычислялись точно. На практике, даже если коэффициенты и
правая часть исходного уравнения имеют простые аналитические выражения, коэффициенты системы уравнений вычисляются с использованием квадратурных формул.
Для построения дискретных задач на последовательности сеток используется следующий подход. Вначале строится система уравнений метода Бубнова—Галерина с использованием кусочно-линейных элементов на треугольниках на самой мелкой сетке.
При этом для вычисления элементов матрицы и вектора правой части применяется
квадратурная формула, точная для многочленов степени ≤ 1. Далее элементы матрицы системы уравнений на более редкой сетке вычисляются как линейные комбинации
элементов матрицы, соответствующей более мелкой сетке. Другими словами, матрицы
систем уравнений на соседних сетках связаны с помощью операторов интерполяции
и проектирования, как и в случае, когда элементы матриц вычисляются точно. Это
свойство дает некоторые удобства при реализации метода на ЭВМ.
В данной работе рассматриваются каскадный алгоритм и полный многосеточный алгоритм на основе несимметричного W -цикла. Доказана оптимальная вычислительная
сложность обоих алгоритмов, т. е. обоснована оценка S ≤ cN , где S — число арифметических операций для достижения порядка точности, совпадающего с порядком
точности решения дискретной задачи с учетом численного интегрирования, N — число
неизвестных системы Бубнова—Галеркина, c — константа, не зависящая от N .
1. Формулировка дифференциальной задачи
В выпуклом открытом многоугольнике Ω ⊂ R2 с границей Γ рассмотрим задачу Дирихле
−
2
P
∂i (aij ∂j u) + au = f в Ω,
(1)
i,j=1
u = 0 на Γ,
(2)
где коэффициенты и правая часть уравнения (1) удовлетворяют условиям гладкости и
эллиптичности:
∂i aij ∈ Lq (Ω), q > 2, i, j = 1, 2; a12 = a21 на Ω̄;
2
2
2
P
P
P
ξi2 на Ω̄ ∀ ξi ∈ R,
aij ξi ξj ≤ ν
µ ξi2 ≤
i=1
i,j=1
ν ≥ µ > 0;
(3)
i=1
a ≥ 0 на Ω̄.
a, f ∈ L2 (Ω);
Тогда задача (1)–(2) имеет единственное решение в W22 (Ω), удовлетворяющее оценке [11]:
kuk2,Ω ≤ c1 |f |0,Ω .
◦
Перейдем к обобщенной формулировке для задачи (1)–(2): найти u ∈ W21 (Ω), такую, что
◦
L(u, v) = f (v) ∀ v ∈ W21 (Ω),
(4)
Многосеточные итерационные алгоритмы в методе конечных элементов...
36
где билинейная форма L и функционал f соответственно определяются соотношениями
Ã
!
2
R P
L(u, v) =
aij ∂j u∂i v + auv dx,
(5)
Ω
i,j=1
f (v) =
Z
f vdx.
(6)
Ω
При условиях (3) задача (4) имеет также единственное решение.
2. Формулировка дискретной задачи
Для построения схемы метода конечных элементов проведем разбиение области Ω̄. Вначале разобьем ее на небольшое количество треугольников так, чтобы получившаяся триангуляция T0 была согласованной, т. е. любые два треугольника должны иметь общую
сторону либо общую вершину или не иметь общих точек. Обозначим через h0 максимальную из длин сторон полученных треугольников. Для i = 1, . . . , l положим Ni = 2i ,
hi = h0 /Ni и разобьем каждый треугольник на Ni2 равных треугольников.
Обозначим
множество всех вершин полученной триангуляции Ti через Ω̄i и введем
T
Ωi = Ω̄i Ω. Через ni обозначим число точек множества Ωi . Каждому узлу y ∈ Ωi
◦
поставим в соответствие базисную функцию ϕiy ∈ W21 (Ω), которая линейна на каждом
треугольнике из Ti , равна единице в узле y и нулю во всех остальных узлах из Ω̄i .
Обозначим через H i линейную оболочку этих функций, т. е. H i = span{ϕiy }, y ∈ Ωi .
◦
Рассматривая (4) на подпространстве H i ⊂ W21 (Ω), получим дискретную задачу:
найти ṽi ∈ H i такую, что
L(ṽi , v) = f (v) ∀ v ∈ H i .
(7)
Обозначим через Mi пространство размерности ni , состоящее из векторов w с компонентами w(x), x ∈ Ωi . Тогда задача (7) эквивалентна системе линейных алгебраических
уравнений
Li v i = f i ,
(8)
где vi ∈ Mi — вектор неизвестных с компонентами vi (y), y ∈ Ωi ; f i ∈ Mi — вектор
правой части с компонентами fi (x) = f (ϕix ), x ∈ Ωi ; Li — матрица размером ni × ni с
элементами Li (x, y) = L(ϕix , ϕiy ), x, y ∈ Ωi .
Произвольному вектору v ∈ Mi поставим в соответствие восполнение в H i :
X
ṽ(x) =
(9)
v(y)ϕiy (x), x ∈ Ω.
y∈Ωi
Из определения базисных функций следует, что
v(y) = ṽ(y),
y ∈ Ωi .
(10)
Таким образом, мы определили изоморфизм между векторами v ∈ Mi и функциями
ṽ ∈ H i .
Введем энергетическую норму для функций
|||v|||Ω = L(v, v)1/2 ,
v ∈ W21 (Ω).
37
Л. В. Гилева
◦
Для функций из W21 (Ω) она эквивалентна норме k · k1,Ω [12]:
◦
c2 kvk1,Ω ≤ |||v|||Ω ≤ c3 kvk1,Ω ,
v ∈ W21 (Ω).
(11)
Введем также скалярное произведение и нормы для векторов:
(v, w)i =
X
x∈Ωi
1/2
1/2
v(x)w(x), kvki = (v, v)i , |||v|||i = (v, Li v)i , v, w ∈ Mi .
На основании (9) и (10) для изоморфной пары v ∈ Mi , ṽ ∈ H i имеем
(12)
|||v|||i = |||ṽ|||Ω .
Норма |ṽ|0,Ω эквивалентна норме kvki с множителем hi [13]:
c4 hi kvki ≤ |ṽ|0,Ω ≤ c5 hi kvki .
(13)
Подпространства H i обладают свойством вложенности H i−1 ⊂ H i , т. е. любая базисная функция ϕi−1
∈ H i−1 является линейной комбинацией базисных функций ϕiy ∈ H i :
x
ϕi−1
=
x
X
βxy ϕiy ,
y∈Ωi
x ∈ Ωi−1 .
(14)
Введем оператор интерполяции Ii−1 : Mi−1 → Mi . Пусть v — произвольный вектор
из Mi−1 . Тогда компоненты вектора w = Ii−1 v ∈ Mi определяются по формуле
w(y) =
X
v(x)βxy ,
x∈Ωi−1
y ∈ Ωi ,
(15)
с коэффициентами βxy из (14). Отметим, что восполнения векторов v и w совпадают,
т. е. ṽ = w̃. Таким образом, оператор Ii−1 соответствует тождественному оператору на
подпространстве H i−1 .
Введем также оператор проектирования Ri : Mi → Mi−1 . Для вектора w ∈ Mi
компоненты вектора v = Ri w ∈ Mi−1 вычисляются по формуле
v(x) =
X
y∈Ωi
βxy w(y),
x ∈ Ωi−1 .
(16)
∗
Из (15) и (16) следует, что (Ri w, v)i−1 = (w, Ii−1 v)i , v ∈ Mi−1 , w ∈ Mi , т. е. Ri = Ii−1
.
При условиях (3) задача (7) имеет единственное решение, которое удовлетворяет
оценкам [12, 13]
|||u − ṽi |||Ω ≤ c6 hi |f |0,Ω ,
|u − ṽi |0,Ω ≤ c7 h2i |f |0,Ω .
(17)
(18)
Многосеточные итерационные алгоритмы в методе конечных элементов...
38
3. Использование численного интегрирования
На практике, даже если функции aij , a и f имеют простые аналитические выражения,
элементы матрицы Li и вектора f i в (8) редко вычисляются точно. Вместо этого они
аппроксимируются с помощью численного интегрирования.
Мы начнем с построения дискретной задачи с учетом численного интегрирования
на самой мелкой сетке. Используя (5) и (6), можем записать
!
Ã
2
R
P
P
(19)
aij ∂j ϕlx ∂i ϕly + aϕlx ϕly dx,
Ll (x, y) = L(ϕlx , ϕly ) =
ω∈Tl ω
fl (x) = f (ϕlx ) =
i,j=1
P R
ω∈Tl ω
f ϕlx dx,
x, y ∈ Ωl .
(20)
Пусть ω — произвольный треугольник из Tl . Квадратурная схема состоит в замене
M
R
P
αm,ω g(bm,ω ), где αm,ω — веса, а точки
интеграла g(x)dx конечной суммой вида
m=1
ω
bm,ω — узлы квадратурной формулы. Будем использовать следующую квадратурную
формулу:
Z
3
meas(ω) X
g(x)dx ≈
g(bm,ω ),
(21)
3
m=1
ω
где meas(ω) — площадь треугольника ω, а bm,ω , m = 1, 2, 3, — вершины ω. Формула
(21) точна для многочленов степени ≤ 1.
Введем аппроксимирующую билинейную форму
L̂(u, v) =
3 X
2
X meas(ω) X
(aij ∂j u∂i v + auv)(bm,ω )
3
m=1 i,j=1
ω∈T
(22)
l
и аппроксимирующий функционал
fˆ(v) =
3
X meas(ω) X
(f v)(bm,ω ).
3
m=1
ω∈T
(23)
l
Подставляя (21) в (19) и (20), получим вместо матрицы Ll матрицу L̂l с элементами
L̂l (x, y) = L̂(ϕlx , ϕly ), x, y ∈ Ωl , а вместо вектора правой части f l вектор f̂ l с компонентами
fˆl (x) = fˆ(ϕlx ), x ∈ Ωl . В итоге мы получаем систему уравнений L̂l wl = f̂ l , которая
эквивалентна следующей задаче: найти w̃l ∈ H l такую, что
L̂(w̃l , v) = fˆ(v) ∀ v ∈ H l ,
где билинейная форма L̂ и функционал fˆ соответственно заданы равенствами (22)
и (23).
Перейдем к построению систем уравнений на более крупных сетках. Матрицы систем
уравнений Бубнова—Галеркина на соседних сетках связаны соотношением [13] Li−1 =
Ri Li Ii−1 , что позволяет по матрице Ll , соответствующей самой мелкой сетке, строить
матрицы Li , i < l, с помощью простых операций суммирования строк и столбцов. Мы
39
Л. В. Гилева
потребуем, чтобы аналогичное свойство выполнялось для матриц L̂i . Предположим, что
на сетке Ωi мы имеем систему уравнений
L̂i wi = f̂ i ,
(24)
где L̂i — матрица размером ni × ni с элементами L̂i (x, y) = L̂(ϕix , ϕiy ), x, y ∈ Ωi , а
f̂ i ∈ Mi — вектор с компонентами fˆi (x) = fˆ(ϕix ), x ∈ Ωi . Положим
L̂i−1 = Ri L̂i Ii−1 ,
f̂ i−1 = Ri f̂ i .
(25)
(26)
Систему уравнений на сетке Ωi−1 определим следующим образом:
L̂i−1 wi−1 = f̂ i−1 .
(27)
Используя равенство (26), определение оператора проектирования (16), линейность fˆ
и соотношение (14), для компонент вектора правой части (27) получим
fˆi−1 (x) = fˆ(ϕi−1
x ), x ∈ Ωi−1 .
(28)
Используя (14), (22), определения операторов интерполяции и проектирования, равенство (25) и учитывая вид элементов матрицы L̂i , получим, что элементы матрицы
L̂i−1 определяются соотношением
i−1
L̂i−1 (x, y) = L̂(ϕi−1
x , ϕy ), x, y ∈ Ωi−1 .
(29)
Из (28) и (29) следует, что для любого i = 0, . . . , l система уравнений (24) эквивалентна
следующей задаче: найти w̃i такую, что
L̂(w̃i , v) = fˆ(v) ∀ v ∈ H i ,
(30)
где билинейная форма L̂ и функционал f˜ соответственно заданы соотношениями (22)
и (23).
Отметим, что задача (30) может быть получена из (7) с использованием составной
квадратурной формулы:
Z
ωi
g(x)dx ≈
3
X meas(ωl ) X
g(bm,ωl ),
3
ω ∈T
m=1
(31)
l
l
ωl ⊂ωi
где ωi — произвольный треугольник триангуляции Ti . Очевидно, что квадратурная схема (31) точна для многочленов степени ≤ 1.
Так как квадратурные формулы (21) и (31) точны для многочленов степени 0, то
решение задачи (30) удовлетворяет оценке [12]
ku − w̃i k1,Ω ≤ c8 hi .
(32)
Кроме того, поскольку квадратурные формулы (21) и (31) точны и для многочленов
степени 1, имеет место оценка в норме L2 [12]
|u − w̃i |0,Ω ≤ c9 h2i .
(33)
Многосеточные итерационные алгоритмы в методе конечных элементов...
40
В квадратурной формуле (21) все веса строго положительны, и она точна для многочленов степени 0. Поэтому существуют константы c10 и c11 , не зависящие от hl , такие,
что [14]
L̂(v, v) ≥ c10 kvk21,Ω ,
(34)
l
L̂(v, w) ≤ c11 kvk1,Ω kwk1,Ω , v, w ∈ H .
(35)
В силу вложенности подпространств H i оценки (34) и (35) справедливы для любых
v, w ∈ H i . Введем норму для функций из H i :
[[v]]Ω = L̂(v, v)1/2 , v ∈ H i .
На основании (34) и (35) она эквивалентна норме k · k1,Ω :
c12 kvk1,Ω ≤ [[v]]Ω ≤ c13 kvk1,Ω , v ∈ H i .
(36)
Отсюда с учетом (11) следует, что
c14 |||v|||Ω ≤ [[v]]Ω ≤ c15 |||v|||Ω , v ∈ H i ,
(37)
т. е. для функций из H i нормы ||| · |||Ω и [[·]]Ω эквивалентны.
На основании (31) матрица L̂i положительно определена, поэтому мы можем ввести
норму для векторов
1/2
[[v]]i = (v, L̂i v)i , v ∈ Mi .
Используя (9) и (10), легко показать, что для изоморфной пары v ∈ Mi , ṽ ∈ H i
(38)
[[v]]i = [[ṽ]]Ω .
Из (37), (12) и (38) следует эквивалентность норм ||| · |||i и [[·]]i , т. е.
c14 |||v|||i ≤ [[v]]i ≤ c15 |||v|||i ,
v ∈ Mi .
(39)
Нам потребуется оценка собственных чисел матрицы L̂i .
Лемма 1. Собственные числа матрицы L̂i удовлетворяют неравенству
0 < λ̂i ≤ c16 ,
i = 0, . . . , l.
(40)
Доказательство. Для собственных чисел матрицы Li мы имеем оценку [13]
0 < λi ≤ c17 ,
i = 0, . . . , l.
(41)
Обозначим через λ∗i и λ̂∗i максимальные собственные числа матриц Li и L̂i соответственно. На основании соотношения Рэлея и оценок (39) и (41) можем записать
неравенства
(v, L̂i v)i
[[v]]2i
λ̂∗i = sup
= sup
≤
(v, v)i
v∈Mi ,
v∈Mi , (v, v)i
v6=0
v6=0
41
Л. В. Гилева
≤ sup
v∈Mi ,
v6=0
c215 |||v|||2i
(v, Li v)i
= c215 λ∗i ≤ c215 c17 ,
= c215 sup
(v, v)i
(v, v)
v∈Mi ,
v6=0
т. е. верхняя оценка в (40) справедлива с константой c16 = c215 c17 . Нижняя оценка следует
из положительной определенности матрицы L̂i .
¤
Итак, на последовательности сеток Ωi мы получили последовательность задач: для
заданного f̂ i ∈ Mi найти вектор wi ∈ Mi такой, что
L̂i wi = f̂ i .
(42)
Наша цель состоит в решении задачи (42) при i = l. Для этого мы используем два
варианта многосеточных методов (каскадный алгоритм и многосеточный алгоритм с
использованием несимметричного W -цикла) и докажем их сходимость.
4. Формулировка каскадного алгоритма
Каскадный алгоритм состоит в последовательном решении задач (42). При i = 0 число
уравнений в (42) невелико и задача решается прямым методом. На более мелких сетках
приближенные решения получают итерационным методом. Сформулируем каскадный
алгоритм с некоторым абстрактным итерационным процессом Si (сглаживателем).
Каскадный алгоритм:
1)
2)
u0 = L̂−1
0 f̂ 0 ;
для i = 1, 2, . . . , l цикл {
2.1. zi = Ii−1 ui−1 ;
2.2. полагаем ui = Si (L̂i , zi , f̂ i ); }.
(43)
В качестве сглаживателя рассмотрим два итерационных процесса.
Метод сопряженных градиентов (mi итераций); процедура Si (L̂i , zi , f̂ i ):
3)
4)
5)
y0 = zi ; p0 = r0 = f̂ i − L̂i y0 ; σ0 = (r0 , r0 )i ;
для i = 1, 2, . . . , mi цикл {
если σk−1 = 0, то ymi = yk−1 и перейти на 5;
αk−1 = σk−1 /(pk−1 , Li pk−1 )i ;
yk = yk−1 + αk−1 pk−1 ;
rk = rk−1 − αk−1 L̂i pk−1 ;
σk = (rk , rk )i ; βk = σk /σk−1 ;
pk = rk + βk pk−1 ; };
полагаем Si = ym .
(44)
Метод простых итераций (mi итераций); процедура Si (L̂i , zi , f̂ i ):
3)
4)
5)
y0 = zi ;
для k = 1, 2, . . . , mi цикл {
π(2k + 1)
;
τk−1 = (Λ∗i )−1 cos−2
2(2mi + 1)
yk = yk−1 − τk−1 (L̂i yk−1 − f̂ i )i };
полагаем Si = ymi .
(45)
Многосеточные итерационные алгоритмы в методе конечных элементов...
42
Здесь Λ∗i — верхняя оценка собственных чисел λ̂ оператора L̂i в пространстве Mi , т. е.
L̂i ϕ = λ̂ϕ. В (45) она нужна в явной форме, поэтому она находится из леммы Гершгорина [15] и удовлетворяет неравенству
λ̂∗i = max λ̂ ≤ Λ∗i ≤ c18 max λ̂ = c18 λ̂∗i .
(46)
λ̂∈Sp (L̂i )
λ̂∈Sp (L̂i )
В (44) в явном виде она не используется, поэтому в теоретических выкладках ее можно
положить равной λ̂∗i , т. е. (46) тогда выполняется с константой c18 = 1.
Зафиксируем целое i ∈ [1, l] и определим линейный оператор Bi : Mi → Mi , ставящий
в соответствие погрешности начального приближения δ 0 = wi − zi для решения задачи
(42) погрешность конечной аппроксимации δ 1 = wi − ui , полученной на i-м уровне
каскадного алгоритма (43) с помощью метода сопряженных градиентов (44), т. е. δ 1 =
Bi δ 0 .
Введем норму для оператора B : Mi → Mi :
[[B]]i = sup
v∈Mi ,
v6=0
[[Bv]]i
.
[[v]]i
Оператор Bi является матричным полиномом от L̂i , т. е.
Bi = I +
mi
X
γk L̂ki .
(47)
k=1
Он обладает следующим важным свойством [16]. При фиксированном начальном приближении y0 (и при фиксированной погрешности начального приближения δ 0 ) метод
сопряженных градиентов минимизирует норму погрешности [[δ 1 ]]i для конечного приближения ymi среди всех многочленов вида (47) с произвольными коэффициентами γk .
Аналогично определим оператор Qi : Mi → Mi для метода простых итераций (45),
т. е. δ 1 = Qi δ 0 . Оператор Qi является матричным полиномом и имеет вид
Qi =
m
i −1
Y
k=0
(I − τk L̂i ).
Для него выполняются следующие оценки [3]:
q
λ̂∗i
kwki ,
[[Qi w]]i ≤
2mi + 1
[[Qi w]]i ≤ [[w]]i ∀ w ∈ Mi .
(48)
(49)
(50)
5. Сходимость каскадного алгоритма
Теорема 1. Пусть оператор L̂i задачи (42) является самосопряженным, положительно определенным и удовлетворяет условию (25). Тогда для приближенного решения ui , полученного на i-м уровне каскадного алгоритма (43) с одним из итерационных
процессов (44) или (45), справедлива оценка
43
Л. В. Гилева
[[wi − ui ]]i ≤ c19
i
X
j=1
hj−1
.
2mj + 1
(51)
Доказательство. Обозначим через ei погрешность на i-м уровне каскадного алгоритма, т. е.
ei = wi − ui , i = 0, 1, . . . , l.
Для метода сопряженных градиентов (44) имеем ei = Bi (wi − zi ), а для метода
простых итераций (45) — ei = Qi (wi − zi ). Оператор Qi имеет вид (47), поэтому на
основании приведенного выше свойства оператора Bi получаем
[[ei ]]i ≤ [[Qi (wi − zi )]]i .
На основании (43) можем записать
[[ei ]]i ≤ [[Qi (wi − Ii−1 wi−1 )]]i + [[Qi Ii−1 (wi−1 − ui−1 )]]i =
= [[Qi (wi − Ii−1 wi−1 )]]i + [[Qi Ii−1 ei−1 ]]i .
(52)
Оценим первое слагаемое в правой части (52). Для этого воспользуемся (49), (13),
(40) и (33):
q
q
∗
λ̂i
λ̂∗i −1
kwi − Ii−1 wi−1 ki ≤ c−1
h |w̃i − w̃i−1 |0,Ω ≤
[[Qi (wi − Ii−1 wi−1 )]]i ≤
4
2mi + 1
2mi + 1 i
q
1/2
λ̂∗i −1
c−1
4 c9 c14
−1
≤ c4
h (|w̃i − u|0,Ω + |u − w̃i−1 |0,Ω ) ≤
h−1 (h2i + h2i−1 ).
2mi + 1 i
2mi + 1 i
Из построения сеток следует, что hi = hi−1 /2. Поэтому
1/2
[[Qi (wi − Ii−1 wi−1 )]]i ≤
5 c−1
4 c9 c14
hi−1 .
2 2mi + 1
(53)
Для оценки второго слагаемого в правой части (52) используем (50):
1/2
[[Qi Ii−1 ei−1 ]]i ≤ [[Ii−1 ei−1 ]]i = (L̂i Ii−1 ei−1 , Ii−1 ei−1 )i
=
1/2
1/2
= (Ri L̂i L̄i−1 ei−1 , ei−1 )i−1 = (L̂i−1 ei−1 , ei−1 )i−1 = [[ei−1 ]]i−1 .
Объединим эту оценку с (52) и (53):
1/2
5 c−1 c9 c14
[[ei ]]i ≤ 4
hi−1 + [[ei−1 ]]i−1 .
2 2mi + 1
(54)
Далее для доказательства используем индукцию по i. При i = 0 имеем [[e0 ]]1 = 0. Из
(54) при i = 1 следует, что
1/2
5 c−1
4 c9 c14
[[e1 ]]i ≤
h0 ,
2 2m1 + 1
5
1/2
т. е. (51) справедливо при i = 1 с константой c19 = c−1
4 c9 c14 . Предположим, что (51)
2
справедливо на уровне i − 1:
[[ei−1 ]]i−1 ≤ c19
i−1
X
j=1
hj−1
.
2mj − 1
Многосеточные итерационные алгоритмы в методе конечных элементов...
44
Объединяя это неравенство с (54), получаем, что оценка (51) выполняется с константой
5
1/2
c19 = c−1
4 c9 c14 .
2
¤
Путем простого анализа последовательности вычислений с учетом разреженности
матриц L̂i устанавливается верхняя оценка числа арифметических операций в каскадном алгоритме:
l
X
Sl ≤ c20
(mj + c21 )nj
(55)
j=1
с константами c20 , c21 , не зависящими от nj и mj , но разными для итерационных процессов (44) и (45). Ясно, что для последнего процесса эти константы существенно меньше.
Число итераций m1 , . . . , ml−1 выберем так, чтобы минимизировать Sl как функцию
от m1 , . . . , ml−1 при фиксированной величине правой части (51) при i = l [3]. В результате mi выбирается как наименьшее целое, удовлетворяющее неравенству
µ
¶1/2
nl hi−1
2mi + 1 ≥ (2ml + 1)
.
(56)
ni hl−1
Теорема 2. Погрешность каскадного алгоритма (43) с mi итерациями метода
сопряженных градиентов (44) или метода простых итераций (45), где mi выбирается
из условия (56) при фиксированном ml , оценивается как
|||ul − wl |||l ≤ c22
hl
.
2ml + 1
Для восполнения ũl ∈ H l вектора ul справедлива оценка
¶
µ
c22
hl .
|||ũl − u|||Ω ≤ c3 c8 +
2ml + 1
(57)
(58)
Число арифметических операций оценивается сверху величиной
Sl ≤ (c23 ml + c24 )nl .
(59)
Доказательство. Для задачи Дирихле имеем оценку
ni ≤ 4i−l nl .
(60)
Из построения сеток следует равенство
hi = 2l−i hl .
(61)
С учетом этим соотношений можно переписать (56) в виде
2mi + 1 ≥ (2ml + 1)23(l−i)/2 .
(62)
Используя это неравенство вместе с (61) в (51) при i = l и оценивая сумму геометрической прогрессии, получим
√
l
l
X
X
hl
hl
2 2
hj−1 3(j−l)/2
(j−l)/2
c19
2
= 2c19
.
2
≤√
[[wl − ul ]]l ≤ c19
2m
+
1
2m
+
1
2m
+
1
2
−
1
l
l
l
j=1
j=1
45
Л. В. Гилева
Привлекая эквивалентность норм (39), мы приходим к оценке (57) с константой
√
2 2 −1
c14 c19 .
c22 = √
2−1
На основании неравенства треугольника, эквивалентности норм (12) и (11) и оценки
(32) мы получаем цепочку неравенств
|||ũl − u|||Ω ≤ |||ũl − w̃l |||Ω + |||w̃l − u|||Ω ≤ |||ul − wl |||l + c3 kw̃l − uk1,Ω ≤ |||ul − wl |||l + c3 c8 hl .
Вместе с уже доказанной оценкой (57) это приводит к (58).
Для оценки числа арифметических операций напомним, что mi выбирается как наименьшее целое, удовлетворяющее (62). Поэтому
1
1
mi ≤ (2ml + 1)23(l−i)/2 + .
2
2
Используя это неравенство и соотношение (60) и оценивая сумму геометрической прогрессии, из (55) получим
l µ
X
1
¶
1
Sl ≤ c20
(2ml + 1)2
+ + c21 4j−l nl =
2
2
j=1
¶X
¶X
µ
µ
l
l
´
1
1
(j−l)/2
j−l
(2ml + 1)
+ c21
2
+
4
nl ≤
= c20
2
2
j=1
j=1
Ã
√
µ
¶!
2
1
4
1
√ +
≤ c20 √
nl ,
c21 +
ml + c20
2
2−1
2− 2 3
3(l−j)/2
т. е. справедлива оценка (59) с константами
√
µ
¶
1
4
1
2
√ +
и c24 = c20
c21 +
.
c23 = c20 √
2
2−1
2− 2 3
¤
Из оценки (58) вытекает, что число итераций на самом верхнем уровне следует выбирать из условия c3 c8 ≈ c22 /(2ml + 1). Тогда погрешность итерационного процесса приобретает такой же порядок малости, как и погрешность дискретной задачи с учетом
численного интегрирования. Хотя константы c3 , c8 и c22 неизвестны, видна независимость ml от числа уровней и количества неизвестных.
6. Формулировка многосеточного алгоритма с использованием
несимметричного W -цикла
Начнем с формулировки алгоритма несимметричного W -цикла для решения задачи (42)
(M Gi -алгоритма). Предположим, что мы имеем некоторое начальное приближение z0i ∈
Mi решения задачи (42). В результате получим приближенное решение z1i = M Gi (z0i , f̂ i ).
Алгоритм M Gi
Многосеточные итерационные алгоритмы в методе конечных элементов...
46
Если i = 0, то полагаем
z10 = M G0 (z0i , f̂ 0 ) = L̂−1
0 f̂ 0
и алгоритм M Gi завершен.
Если i > 0, то алгоритм состоит из нескольких шагов.
1. Положим u∗0 = z0i и вычислим
u∗k = u∗k−1 − τk−1 (L̂i u∗k−1 − f̂ i ),
где
τk−1 =
k = 1, . . . , m,
1 + cos α
,
− cos(2k + 1)α)
α=
Λ∗i (cos α
π
.
2m + 2
(63)
2. Положим gi−1 = Ri (L̂i u∗m − f̂ i ), где Ri : Mi → Mi−1 — оператор проектирования.
3. Возьмем z∗0 = 0 ∈ Mi−1 и повторим алгоритм M Gi−1 2 раза:
z∗S = M Gi−1 (z∗s−1 , gi−1 ),
s = 1, 2.
4. Положим
z1i = M Gi (z0i , f̂ i ) = u∗m − Ii−1 z∗2 .
Теперь сформулируем полный многосеточный алгоритм для решения последовательности задач (42).
Алгоритм F M G
1)
2)
u0 = L̂−1
0 f̂ 0 ;
для i = 1, 2, . . . , l цикл {
2.1. ui = Ii−1 ui−1 ;
2.2. повторить алгоритм M Gi t раз :
ui := M Gi (ui , f̂ i ).
В результате мы получим ul — приближенное решение задачи (42) при i = l.
Рассмотрим матричный полином (48), где параметры τk заданы соотношениями (63).
Для него выполняется оценка [13]
kL̂i Qi ki ≤ λ̂∗i
tg α/2
,
m+1
π
.
2m + 2
α=
(64)
7. Сходимость алгоритма F M G
Лемма 2. Существует константа c25 такая, что имеет место оценка
c25
−1
kL̂−1
i − Ii−1 L̂i−1 Ri ki ≤
.
(65)
λ̂∗i
Доказательство. Пусть gi ∈ Mi — произвольный вектор. Рассмотрим систему
уравнений
Li v i = g i .
(66)
Любому вектору gi ∈ Mi соответствует функция g ∈ H i такая, что (66) является системой Бубнова—Галеркина для задачи [13]:
◦
L(v, w) = (g, w)Ω ∀ w ∈ W21 (Ω),
47
Л. В. Гилева
т. е. система (66) эквивалентна задаче
L(ṽi , w) = (g, w)Ω ∀ w ∈ H i .
При этом выполняется оценка (18), т. е.
|v − ṽi |0,Ω ≤ c7 h2i |g|0,Ω .
(67)
Рассмотрим также систему уравнений
L̂i w∗i = gi ,
которая эквивалентна задаче
L̂(w̃i∗ , w) = (g, w)Ω ∀ w ∈ H i .
Очевидно, что для w̃i∗ справедлива оценка вида (33), т. е.
|v − w̃i∗ |0,Ω ≤ c9 h2i .
(68)
С помощью неравенства треугольника и оценок (67) и (68) получаем, что
|ṽi − w̃i∗ |0,Ω ≤ |ṽi − v|0,Ω + |v − w̃i∗ |0,Ω ≤ (c7 |g|0,Ω + c9 )h2i ≤ c26 h2i .
С учетом эквивалентности норм (13) для изоморфных векторов получим оценку
kvi − w∗i ki ≤ c4−1 c26 hi .
(69)
Используя неравенство треугольника, можем записать неравенство
−1
−1
−1
kL̂−1
i − Ii−1 L̂i−1 Ri ki ≤ kL̂i − Li ki +
−1
−1
−1
+kL−1
i − Ii−1 Li−1 Ri ki + kIi−1 (Li−1 − L̂i−1 )Ri ki .
(70)
Оценим первое слагаемое в правой части (70). На основании определения матричной
нормы имеем
−1
k(L̂−1
i − Li )gi ki
−1
=
kL̂−1
−
L
k
=
sup
.
i
i
i
kgi ki
gi ∈Mi ,
gi 6=0
Пусть ḡi ∈ Mi — вектор, на котором достигается супремум в правой части последнего
равенства. Ему соответствуют векторы v̄i и w̄∗i , являющиеся решениями систем уравнений Li v̄i = ḡi и L̂i w̄∗i = ḡi соответственно. Тогда, используя оценку (69), получим
неравенства
kL̂−1
i
−
L−1
i ki
−1
k(L̂−1
kv̄i − w̄∗i ki
c−1
c26 hi
i − Li )ḡi ki
=
=
≤ 4
.
kḡi ki
kḡi ki
kḡi ki
Отсюда следует, что при фиксированном ḡi существует константа c27 такая, что
−1
kL̂−1
i − Li ki ≤ c27 .
(71)
Многосеточные итерационные алгоритмы в методе конечных элементов...
48
Наконец, с помощью (40) получим, что
−1
kL̂−1
i − Li ki ≤
c16 c27
λ̂∗i
(72)
.
Для второго слагаемого в правой части (70) имеем оценку [13]
−1
kL−1
i − Ii−1 Li−1 Ri ki ≤
c28
.
λ∗i
При доказательстве леммы 1 мы, в частности, получили оценку λ̂∗i ≤ c215 λ∗i . Используя
ее, можем записать:
c215 c28
−1
−1
kLi − Ii−1 Li−1 Ri ki ≤
.
(73)
λ̂∗i
Нам осталось оценить третье слагаемое в правой части (70). Пусть v ∈ Mi−1 — произвольный вектор. На основании эквивалентности норм (13) можем получить неравенство
kIi−1 vki
≤ 2c−1
4 c5 .
kvki−1
(74)
Норма оператора интерполяции определяется соотношением
kIi−1 kMi ←Mi−1 = sup
v∈Mi−1 ,
v6=0
kIi−1 vki
.
kvki−1
С учетом (74) получаем оценку
kIi−1 kMi ←Mi−1 ≤ 2c−1
4 c5 .
(75)
Из определения оператора проектирования следует равенство
(Ii−1 v, w)i = (v, Ri w)i−1 ∀ v ∈ Mi−1 ,
w ∈ Mi .
(76)
Положим v = Ri w. Тогда с помощью неравенства Коши и оценки (75) для левой части
(76) получим оценку
|(Ii−1 v, w)i | ≤ kIi−1 vki kwki ≤ kIi−1 kMi ←Mi−1 kvki−1 kwki ≤ 2c−1
4 c5 kRi wki−1 kwki .
(77)
В результате из (76), (77) имеем
kRi wki−1
≤ 2c−1
4 c5 .
kwki
(78)
На основании определения нормы оператора проектирования
kRi kMi−1 ←Mi = sup
w∈Mi ,
w6=0
kRi wki−1
kwki
из (78) в силу произвольности w следует оценка
kRi kMi−1 ←Mi ≤ 2c−1
4 c5 .
(79)
49
Л. В. Гилева
Используя оценку (71), которая справедлива для любого i, и оценку (40), получим
−1
kL̂−1
i−1 − Li−1 ki−1 ≤
c16 c27
λ̂∗i
(80)
.
В результате с помощью (75), (79) и (80) получаем оценку
−1
−1
−1
kIi−1 (L̂−1
i−1 − Li−1 )Ri ki ≤ kIi−1 kMi ←Mi−1 kL̂i−1 − Li−1 ki−1 kRi kMi−1 ←Mi ≤
2
4c−2
4 c5 c16 c27
λ̂∗i
. (81)
Наконец, привлекая (72), (73) и (81) в (70), мы приходим к оценке (65) с константой
2
c25 = c16 c27 + c215 c28 + 4c−2
4 c5 c16 c27 .
¤
Объединяя (64) и (65), получим, что выполняется критерий сходимости
−1
kL̂−1
i − Ii−1 L̂i−1 Ri ki kL̂i Qi ki ≤ η(m),
где
η(m) =
c25 π
c25 tg α/2
≈
.
m+1
4(m + 1)2
(82)
(83)
Зафиксируем целое i ∈ [1, l] и определим линейный оператор Bi : Mi → Mi , ставящий
в соответствие погрешности начального приближения δ 0 = wi − z0i для решения задачи
(42) погрешность конечной аппроксимации δ 1 = wi − z1i , полученной алгоритмом M Gi ,
т. е. δ 1 = Bi δ 0 .
Теперь мы можем сформулировать результат, характеризующий сходимость алгоритма M Gi [13].
Теорема 3. Пусть операторы L̂i являются самосопряженными и положительно
определенными. Пусть также выполняется критерий сходимости (82), где функция
η(m) не зависит от i и стремится к нулю при m → ∞. Тогда для любого ξ ∈ (0, 1)
существует m0 такое, что для всех m ≥ m0 :
kBi ki ≤ ξ, i = 1, . . . .l.
(84)
Из оценки (84) вытекает, что для любого ε ∈ (0, 1) существует m0 , не зависящее от hi ,
такое, что для всех m ≥ m0 погрешность приближенного решения z1i задачи (42), полученного алгоритмом M Gi , уменьшается с множителем ε по сравнению с погрешностью
начального приближения z01 , т. е.
kwi − z1i ki ≤ εkwi − z0i ki .
(85)
Отметим, что на основании (83) множитель ε в (85) зависит от числа итераций m следующим образом [13]:
ε = O((m + 1)−2 ).
Рассмотрим реализацию алгоритма F M G. Выберем m так, чтобы выполнялось неравенство ε < 1 для ε из (85). Тогда мы можем выбрать целое t, такое что
εt < c4 c−1
5 /4.
(86)
Приведем результат, характеризующий сходимость алгоритма F M G [13, теорема 4.15].
Многосеточные итерационные алгоритмы в методе конечных элементов...
50
Теорема 4. Пусть m выбрано так, что (85) выполняется с ε < 1 и пусть t удовлетворяет условию (86). Тогда при использовании алгоритма F M G мы получим последовательность приближенных решений ui , i = 1, . . . , l, задач (42), которые удовлетворяют оценке
1+ρ
kwi − ui ki ≤ c9 hi
, ρ = 4εt c−1
4 c5 < 1,
c5 (1 − ρ)
а для их восполнений ũi ∈ H i справедлива оценка
|u − ũi |0,Ω ≤ c9 h2i · 2/(1 − ρ).
Наконец, сформулируем результат, содержащий оценку числа арифметических операций в алгоритме F M G (см. [13, лемма 5.1.4]).
Теорема 5. В алгоритме F M G для достижения точности, обусловленной погрешностью дискретной задачи с учетом численного интегрирования, т. е. для выполнения оценок
|w̃l − ũl |0,Ω ≤ c9 h2l ,
|u − ũl |0,Ω ≤ 2c9 h2l
требуется O(nl ) арифметических операций.
Список литературы
[1] Ильин В.П., Свешников В.М. О разностных методах на последовательности сеток //
Численные методы механики сплошной среды: Информ. бюл. / АН СССР. Сиб. отд-ние.
Отд-ние механики и процессов управления. 1971. Т. 2, № 1. С. 43–55.
[2] Deuflhard P. Cascadic conjugate-gradient methods for elliptic partial differential equations.
I. Algorithm and numerical results // Technical Report SC 93-23. Berlin: Konrad-ZuseZentrum (ZIB), 1993.
[3] Shaidurov V.V. Some estimates of the rate of convergence for the cascadic conjugategradient method // Comp. Math. Applic. 1996. Vol. 31, N 4–5. P. 161–171.
[4] Shaidurov V.V. Some estimates of the rate of convergence for the cascadic conjugategradient method. Magdeburg, 1994 (Preprint Otto-von-Guericke Universität. № 4).
[5] Bornemann F.A., Deuflhard P. The cascadic multigrid method for elliptic problems //
Numer. Math. 1996. Vol. 75. P. 135–152.
[6] Shaidurov V.V., Timmermann G. A cascadic multigrid algorithm for semilinear indefinite
elliptic problems // Computing. 2000. Vol. 64. P. 349–366.
[7] Shaidurov V.V. Cascadic algorithm with nested subspaces in domains with curvilinear
boundary // Advanced Mathematics: Computations and Applications. Novosibirsk:
Computing Center Publisher, 1995. P. 588–595.
[8] Shaidurov V.V., Tobiska L. The convergence of the cascadic conjugate-gradient method
applied to elliptic problems in domain with re-entrant corners // Math. Comput. 2000. Vol. 69,
N 230. P. 501–520.
[9] Гилева Л.В. Каскадный многосеточный алгоритм в методе конечных элементов для
трехмерной задачи Дирихле // Сиб. журн. вычисл. математики. 1998. Т. 1, № 3.
C. 217–226.
51
Л. В. Гилева
[10] Гилева Л.В., Шайдуров В.В. Каскадный многосеточный алгоритм в методе конечных
элементов для трехмерной задачи Дирихле в области с криволинейной границей // Сиб.
журн. вычисл. математики. 2002. Т. 5, № 2. C. 127–147.
[11] Ладыженская О.А., Уральцева Н.Н. Линейные и квазилинейные уравнения эллиптического типа. М.: Наука, 1973. 578 c.
[12] Сьярле Ф. Метод конечных элементов для эллиптических задач: пер. с англ. М.: Мир,
1980.
[13] Shaidurov V.V. Multigrid Methods for Finite Elements. Dordrecht: Kluwer Academic
Publishers, 1995.
[14] Даутов Р.З., Карчевский М.М. Введение в теорию метода конечных элементов: учеб.
пособие. Казань: Изд-во Каз. ГУ, 2004.
[15] Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984. 320 с.
[16] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука,
1978. 592 c.
Поступила в редакцию 11 ноября 2007 г.,
в переработанном виде — 13 мая 2008 г.
Документ
Категория
Без категории
Просмотров
7
Размер файла
253 Кб
Теги
конечный, многосеточные, алгоритм, метод, элементов, итерационные, интегрированный, учетом, численного
1/--страниц
Пожаловаться на содержимое документа