Table Of ContentИнститут проблем механики
Российской Академии Наук
С. Д. Алгазин, Г. И. Кершуков
ЧИСЛЕННЫЕ АЛГОРИТМЫ КЛАССИЧЕСКОЙ
МАТФИЗИКИ.
VII. Об уравнении Пуассона в цилиндре.
Препринт № 722
Москва 2003 г.
Аннотация.
Рассматривается трёхмерное уравнение Пуассона в цилиндре с
неоднородными краевыми условиями и правой частью
обеспечивающими гладкость решения. Для приближенного
нахождения этого решения построен численный алгоритм без
насыщения. Указан эффективный способ решения соответствующей
дискретной задачи.
The summary.
The three-dimensional Poisson equation in the cylinder with non-
uniform regional conditions and with the right of part ensuring smooth of
the solution is considered. For approximate finding of this solution the
numerical algorithm without saturation is constructed. The effective way
of the solution of the appropriate discrete problem is specified.
055(02)2 Институт проблем механики РАН 2003
2
Введение.
Это седьмой препринт серии объявленной в [1]. Предыдущие
препринты этой серии [1-6]. Препринт посвящён описанию
программного комплекса решения задачи об уравнении Пуассона в
цилиндре. Подробное описание алгоритма приведено в [7].
В [7] описана дискретизация уравнения Пуассона в цилиндре с
однородными краевыми условиями. В этой методике требование
однородности краевых условий существенно, поскольку идея
дискретизации состоит в построении конечномерной задачи со
спектральными свойствами, аналогичными дифференциальной
задаче. Вместе с тем практические задачи часто приводят к
уравнению Пуассона с неоднородными краевыми условиями. В
качестве примера можно привести задачу о деформации горного
массива над нефтяным пластом, вызванную падением давления в
пласте после откачки нефти. Если считать горный массив упругим
телом, то объемное расширение удовлетворяет уравнению Пуассона.
В качестве области, в которой это уравнение рассматривается,
можно выбрать цилиндр с достаточно большим радиусом основания,
снося на его поверхность граничные условия из бесконечности.
Ниже построен численный алгоритм без насыщения [7] для
уравнения Пуассона в цилиндре с неоднородными краевыми
условиями. Особенности численных алгоритмов без насыщения
является то, что они в отличие от разностных методов приводят к
дискретной задаче с полностью заполненной матрицей.
Эффективное решение которой и представляет основную трудность.
Ниже подробно описан прием, позволяющий свести обращение
3
матрицы большого размера к обращению нескольких матриц
меньшего размера. Можно сказать, что, вместо решения трехмерной
задачи, решается несколько двухмерных задач. В случае кругового
цилиндра задача еще упрощается и сводится к решению нескольких
одномерных задач.
Описание методики проведём для случая кругового цилиндра. По
ходу изложения будет сказано, как обобщить эти результаты на
случай цилиндра с произвольным основанием.
I. Постановка задачи и дискретизация.
Рассматривается уравнение Пуассона с краевым условием Дирихле
∆u+f=0, (I.1)
u| =g. (I.2)
∂G
Решение ищется в круговом цилиндре единичного радиуса, т.е. в
области G={r≤1;a≤ z ≤b}, где (r,ϑ,z) - цилиндрические координаты.
Для примера рассмотрено краевое условие Дирихле. Проведённые
ниже рассуждения без труда обобщаются на другие типы краевых
условий.
Будем обозначать: g =g (r,ϑ) - граничное условие на донышке
a a
цилиндра (z=a); g =g (r,ϑ) - граничное условие на крышке
b b
цилиндра (z=b); g =g (ϑ, z) - граничное условие на боковой
c c
поверхности цилиндра.
∂2u
∆u=∆ u+ , (I.3)
r,ϑ ∂ z2
где ∆ - плоский оператор Лапласа. Из (I.3) следует
r,ϑ
∂2u 2π
u(ζ,z)= ∫K(ζ,ξ)(f (ξ,z)+ )dξ+ ∫ K (ζ,ϑ)g (ϑ,z)dϑ. (I.4)
- ∂ z2 0 c
|ζ|≤1 0
Здесь
4
1
K(ζ,ξ) = − ln|(1-ζξ)/(ζ−ξ)|, ζ=ρexp(iϕ), ξ=rexp(iϕ)
2π
- функция Грина оператора Лапласа с краевым условием Дирихле
(однородным);
1 1−ρ2
K (ζ,ϑ) = , ζ=ρexp(iϕ)
0 2π1+ρ2 −2ρcos(ϑ−ϕ)
- ядро Пуассона.
Выберем в круге (сечении цилиндра плоскостью z=const) сетку
из m окружностей и N=2n+1 точек на каждой окружности. На
границе круга также выберем N=2n+1 точек. Причём радиус ν-ой
окружности r =cos(2ν−1)/π/4m), ν=1,2,...,m. По окружностям узлы
ν
располагаются через равные углы ϑ = 2πl/N, l=0,1,...,2n.
l
Применим для функции
∂2u
F(ξ,z) =f(ξ,z)+
∂ z2
интерполяционную формулу К. И. Бабенко для функции двух
переменных в круге [8, стр. 212], а для функции g (ϑ, z) применим
c
интерполяцию тригонометрическим многочленом (здесь z -
фиксировано), тогда из (I.4) получаем приближённое значение
~
u(ζ,z) для функции u(ζ,z)
~ 2n
u(ζ,z)=∑H (ζ)F(ξ,z)+∑H0(ζ)g (ϑ,z), (I.5)
i i j c j
i j=0
где
H (ζ) = ∫K(ζ,ξ)l (ξ)dξ,
i i
|ζ|≤1
2 n
H0= (0.5+∑ρl cosl(ϕ−ϑ )),ζ=ρexp(iϕ).
j N j
l=1
5
Здесь l (ξ), i=1,2,...,R, R=mN - фундаментальные функции
i
интерполяционной формулы К. И. Бабенко в круге. Конкретный вид
этих функций приведен в [8, стр. 212]. Для дальнейших рассуждений
он неважен.
Проведём теперь дискретизацию ∂2v/∂z по z. Выберем для v(ξ,z)
интерполяционную формулу по k узлам
~ (−1)k+1(x−1)
v(ξ,z)= T (x)g (ξ)
2 k a
2(1−x2) k 1 k−1 ~
+ ∑ ∑'cos(qψ)T (x)v(ξ,z )
k sin2ψ j q j
j=1 j q=0
(x+1)
+ T (x)g (ξ). (I.6)
2 k b
~
Здесь v(ζ,z) - приближённое значение функции v(ζ,z),
T (x)=cos(karccosx) - многочлен Чебышева степени k; z=((b-
k
a)x+a+b)/2; ψ =(2j−1)π/2k, x =cosψ; z =((b−a)x +a+b)/2,
j j j j j
j=1,2,...,k; штрих у знака суммы означает, что слагаемое при q=0
берётся с коэффициентом ½ .
Обозначим
(−1)k+1(x−1)
L (x)= T (x),
a 2 k
2(1−x2) 1 k−1
L (x)= ∑'co (sqψ )T (x),
j k sin2ψ j q
j q=0
(x+1)
L (x)= T (x).
b 2 k
Дифференцируя (I.6) два раза по z, получим
~
∂2v 4 k ~
= (L'' (x)g (ξ)+∑L'' (x)v(ξ,z )+L'' (x)g (ξ)). (I.7)
∂ z2 (b−a)2 a a j j b b
j=1
6
Обозначим
4
Φ(ζ,z)=∑H (ζ)[f(ζ ,z)+ {L'' (x)g (ξ )
p p (b−a)2 a a p
p (I.8)
2n
+L'' (x)g (ξ )]}+∑H0(ζ)g (ϑ,z),
b b p q c q
q=0
4
D = L'' (x ); i,j = 1,2,...,k
ij (b−a)2 j i
матрицу размера k×k. Тогда из (7.5), (7.7) получаем
k ~
u(ζ ,z )=∑H ∑D u(ζ ,z )+Φ(ζ ,z ). (I.9)
q i q p ij p j q j
p j=1
Здесь H = H (ζ ), где p,q=1,2,...,R пробегают узлы интерполяции в
qp p q
круге, т.е. H - матрица дискретной задачи Дирихле для плоского
оператора Лапласа ; z ,i =1,2,...,k пробегает узлы интерполяции по z.
i
Эффективный способ вычисления элементов матрицы H указан в
[7].
Исследование структуры конечномерной задачи.
Перенумеруем узлы в цилиндре сначала по z, а затем по ζ, т.е.
быстрее всего меняется индекс i, потом q. Тогда имеем из
соотношений (I.9) в матричной форме
u = H ⊗Du+Φ, (I.10)
где u - вектор столбец, содержащий приближённые значения
искомого решения в узлах сетки; Φ - вектор столбец, содержащий
значения соответствующей функции (см. I.8) в узлах сетки; знаком
⊗ обозначено кронекеровское произведение матриц H и D.
Размерность этих векторов N = Rk равна числу внутренних узлов
g
сетки. Решив конечномерную задачу (I.10) получим приближённое
значение решения в узлах сетки. В остальных точках цилиндра
7
решение может быть восстановлено по используемым выше
интерполяционным формулам.
Таким образом, для решения системы линейных уравнений
(I.10) требуется обратить матрицу I − H⊗D большого размера
N xN . Например, для сетки m=5, N=11, k=10 эта матрица 550х550,
g g
и она полностью заполнена. Современные персональные ЭВМ
позволяют работать с такими матрицами, но время обращения
достаточно большое. На ПЭВМ типа Pentium с тактовой частотой
90Мгц время обращения матрицы 550х550 по стандартной
программе занимает примерно 7 минут. Кроме этого, требуется
много памяти для хранения полностью заполненной матрицы. В
этом параграфе описывается как преодолеть эти трудности.
Пусть
k
D=∑λb ,b2 =b , b b =0, q ≠ p (I.11)
q q q q q p
q=1
есть спектральное разложение матрицы D. Такое разложение всегда
можно построить, если D - матрица простой структуры, т.е. имеет
полную систему собственных векторов. Поскольку по построению D
- это матрица дискретного оператора, соответствующего
дифференциальному оператору d2/dz2 с однородными краевыми
условиями при z=a и z=b, то это условие выполняется. Далее в (I.11)
b ,q=1,2,...,k - собственные проекторы на одномерное инвариантное
q
подпространство, λ - соответствующее собственное значение. В
q
практических расчётах размер матрицы D невелик (k=6÷20) и
спектральное разложение можно построить, решив полную
проблему собственных значений для матриц D и D'.
8
k
Заметим, что I = I ⊗I ;I = ∑b , тогда
R k k k
q=1
k k k
I − H⊗D = I ⊗(∑b )− H⊗(∑λb ) = ∑(I −λ H)⊗b .
R q q q R q q
q=1 q=1 q=1
В таком случае имеем следующую формулу для обратной к матрице
I − H⊗D:
k
(I − H⊗D)−1 = ∑(I −λ H)−1 ⊗b , (I.12)
R q q
q=1
которая проверяется непосредственным умножением.
Таким образом, вместо обращения матрицы размера N xN
g g
требуется обратить k матриц размера R×R, т.е. размера, равного
числу узлов сетки в круге. Отметим, что I −λ H - h-матрица и её
R q
обращение сводится к обращению (n+1) матриц размера m×m. Кроме
того, для хранения h-матрицы требуется хранить всего m2(n+1)
элементов. Теперь для того, чтобы решить систему линейных
уравнений (I.10), нужно умножить правую часть соотношения (I.12)
на вектор. Если N = 3µ, µ=1,2,..., то для умножения h-матрица на
вектор можно, для уменьшения количества операций, использовать
быстрое преобразование Фурье. Следовательно, для кругового
цилиндра формула (I.12) даёт исчерпывающее решение
поставленной задачи.
Обсудим теперь изменения, которые следует внести в методику,
чтобы рассмотреть цилиндр с произвольным основанием в виде
области G, с гладкой границей ∂G. Пусть ϕ(z), |z|≤1 конформное
отображение круга на область G. Обозначим
Z =diag(|ϕ'(ζ)|2,...,|ϕ'(ζ )|2). Тогда аналогично получаем, что в
1 R
формулы этого параграфа вместо H нужно подставить матрицу HZ.
Изменится также правая часть Φ соотношения (I.11). Вместо f(ζ ,z)
p
9
в (I.8) войдёт |ϕ'(ζ )|2 f(ζ ,z). Эти изменения приведут к тому, что в
p p
формуле (I.12) придётся численно обращать k матриц размера R×R.
Полученная после обращения матрица будет матрицей общего вида
и, для её умножения на вектор не удастся использовать быстрое
преобразование Фурье.
Обсуждение методики и численный пример.
Для того, чтобы оценить погрешность предложенной методики,
нужно выписать применяемые интерполяционные формулы с
остатком и оценить погрешность дискретизации. Это делается
стандартными приёмами и поэтому не рассматривается.
Исследование устойчивости проводится ниже при рассмотрении
конкретных расчётов. Отметим качественную характеристику
методики. Выше применялась интерполяция решения
многочленами: алгебраическими и тригонометрическими. Известно
[8,9], что интерполяционный многочлен приближает гладкую
функцию тем точнее, чем большим условиям гладкости она
удовлетворяет, т.е. мы получили метод без насыщения. Он
автоматически настраивается на гладкость решения. Гладкость
решения определяется входными данными: граничными условиями
и правой частью.
В качестве численного примера рассматривалось решение
уравнения Пуассона в круговом цилиндре единичного радиуса при
a =0, b=π.
Пусть u(r,ϑ,z) = r3z6cos5ϑ, тогда это решение удовлетворяет
соотношениям (I.1), (I.2), если
10
