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

УДК 517.54:681.3.06

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

В ДВУСВЯЗНЫХ ОБЛАСТЯХ

© 2007 г. В.В. Соболев

Numerical method for stable solution of the mixed boundary problem for harmonic functions in the bounded doubly-connected domains by non classical discrete propounding was elaborated. This method based on preliminary building of conformal mapping of given domain into circle ring and subsequent solution transpose boundary problem for ring by direct methods with procedure regulari-zations. Testing of method with employment of computer programs demonstrates sufficiently high efficiency and precision of it.

1. Постановка задачи

Пусть D - ограниченная двусвязная область с жордановыми кусочно-гладкими связными граничными компонентами Г+ (внешняя) и Г- (внутренняя). На Г+ заданы п+, n+eN пар точек A1+, В+,..., Л+, B ++, расположенных в положительном направ-

п п

+ п+

лении обхода Г . Обозначим Г+= и (Л+ В+),

1 к=1 к' к

Г+ = и (В+, Л+ ,) (Л++ = Л+). Аналогично на Г- зада-

2 к=1 к * к+1' 4 и +1 1 у

ны точки Л- В- ... Л" В- , и"еК, расположенные

1 ' 1 ' ' и ' п

в положительном направлении обхода Г-. Обозначим

„- п -

г- = и(Лк, в,), Г- = и Вк, Лк+1) (Ли-+1 = Л-). Пусть

на Г+ и Ц" дискретно, в отдельных точках е Г+ (у = 1,...,М+), е Г- (у = 1,...,М-) заданы значения действительных непрерывных функций р+ (£), р (О, а на Г^ и Г2- в отдельных точках ^+2) е Г+ (у = 1,...,М2+), е Г22 (У = 1,...,М-) заданы значения действительных непрерывных функций (£), (£): ((£]±1;) = ((, <МС<±2)) = • Требуется найти гармоническую в области D функцию удовлетворяющую граничным условиям:

зн(С(±2})

hZT) = (j = 1.-.Mf),

■Y (j = 1,...,m±).

щегося конформным инвариантом) и равна q = 1/m(D). Пусть % и у± - образы Г± и Г± соответственно при этом отображении, а а11» и a(jt2) - образы точек zf и z<±2): af1) = f (Zf4)' aifi) = f (ZZ±2)). Обозначим arg af1} = afI}, arg af2> = af2).

Пусть z = g (w) - функция, обратная к w = f (z). Рассмотрим композицию функций h o g (w) = H(w). Для нахождения значения h(z) в точке z e D достаточно найти значение H(w) в соответствующей точке w = f (z) e K .

Как известно, при конформном отображении свойство гармоничности функции сохраняется. Поэтому H(w) - гармоническая функция в кольце K. Функция H(w) принимает на у± значения cp(g(a)). На осталь-

+ д ных частях у2 границы кольца имеем _н (a) = Y(a),

дп

где Y(a) = ¥(Z)g\a)\ = y(Z)/\f '(Z)|, Z = g(a) [2].

Таким образом, приходим к следующей смешанной краевой задаче для уравнения Лапласа в кольце K: АН = 0,

Н (af) = ф± (j = 1,..., M± ), (1)

д

дп

H(®(±2)) = Y (j = i,-,m±).

(2)

где

Y *=Y /| / '(^Ч

(3)

дп У

Здесь д / дп обозначает производную в направлении внешней нормали к Г.

2. Редукция к задаче в круговом кольце

Совершим конформное отображение V = /(¿) области D на круговое кольцо к = : д < V < 1} с некоторым д, 0< д < 1, нормированное условием у(2 ) = 1, где 70, еГ+ - произвольно фиксированная точка. Известно [1], что такая функция ^ = у(z) существует и притом единственная. Величина д определяется величиной модуля т(Р) двусвязной области D (являю-

3. Решение редуцированной задачи в круговом кольце методом регуляризации

Будем искать И(л>) в виде Н = Яе Ргде Р(м) - аналитическая в кольце К функция. В общем случае Р^) представляется суммой однозначной аналитической в кольце функции, могущей быть разложенной в ряд Лорана, и логарифма V с действительным коэффициентом:

да

Р(V) = Л • 1п V + ^ (ск - 1ёк)Vк .

Здесь Л С1, ^ (к = 0, ±1, ±2,...) - действительные коэффициенты, которые нужно определить. Полагая

|w| = r, arg w = S, и ограничиваясь достаточно большим натуральным числом N, получаем приближённое равенство (d0 = 0):

N

H(r ■ eiS) s HN (r ■ eiS) = A ■ In r + £ rk (ck cos kS+dk sin kS) =

k=-N

N

= A ■ In r + c0 + £( (rkck + r ~kc-k )cos kS + k=i

+ (rkdk - r-kd_k ) sin kS) , q < r < 1. (4) Граничные условия (1) редуцируются к виду

N

co + £(( + c-k )cos ka]+l) + (dk - d-k )sin ka'f) s pp k=i

(j = 1,..., M+ ), (5)

N

A ■ In q + c0 + £( (qkck + q ~kc-k )cos kaа_1) +

k=i

+ (qkdk - q-kd-k) sin ka(¡X)) s p- (j = MГ ) •

На внешней граничной окружности кольца направление внешней нормали совпадает с направлением по радиусу окружности от её центра, значит,

dH = dH , в то время как на внутренней граничной

dn dr

dH dH

окружности = . Поэтому граничные условия

dn dr (2) редуцируются к виду

N

A + £ k (( - c-k )cos kaf2) + (dk + d-k) sin kaf2)) s *

Записывая необходимые условия минимума

dM

(j = 1,..., M 2+ ),

(6)

N

A / q + £ к ((qkck - ке_к )cos ка(Г2) + к=1

+ (qkdk + q-kd_к)sinка(Г2))/q = (j = 1,-,M-) .

Неизвестные коэффициенты

А, с0, с±1, d±1, ..., c±N, d±N найдём методом регу-ляризующего функционала: мÄ=a + ÄQ^ min . Здесь

а = а(А, Со, С1, С-1, dx, d-1,..., cN, C-N, dN, d-N) =

1 2 = £ ( (exp(ia<+1)) -pp) + j=1

m+ ч2

+ £ ((exp(/aa+2)))/dr- j j=1

M- 2

+£ (((q ■ exp(ia^"1)))-p- ) +

j=1

M-

"£(( (q ■ exp(«Or2)))/dr + j

j=1

— = 0, cA 5c,

дмл= 0, aM,= 0, aM¿ = 0

dck

dd,

суммарная квадратичная невязка уравнений (5), (6); в

качестве «стабилизатора» взят функционал

1 / \ 0 = -1 (И2м(е") + ((ы(еи)/<И))<И +

1 / \ +— | (И2„(д • е") + (йИ м(д • е") / Л) ) Л;

А, Л > 0, - параметр регуляризации.

^к ^"к (к = ± 1, ±2,..., ± N) и учитывая, что

° = 0.(А, с0, с1; с_1; ^ й^,..., ^, с_ ^, dN, й- ^) =

= 2 (А21п2 ц + 2А • с0 • 1п ц + 2с02) +

N Г 2 2

+£ (1 + 42) ( + С-4) +( • + д-ка_к ) +

4=1 ^

+(й4 - к )2+(дк • йк- д- кй- к )2

приходим к следующей системе линейных алгебраических уравнений (СЛАУ):

(верхний значок (0) у всех сумм в формулах (7) означает, что индекс суммирования т не принимает значения 0). Здесь

g¡¡_1 (Л) = М+ + М2+ + М- • 1п2 д + 2Я • 1п д,

g00o(Я) = М_ • 1пд + 2А • 1пд,

g0*1 (л) = М_ • 1п д + 2л • 1п д, g0o* (Л) = М+ + М_ • 1п д + 4Л,

gкm (Л) = gкm + Л • ¿кт (1 + к 2)(1 + д ),

gк-т (Л) = g4-т + 2Л-Якт (1 + к2),

g-кm (Л) = g-кm + 2Л (1 + к'), g-к-т (Л) = g-к-т +Л-5кт (1 + к 2)(1 + д-2к),

g;я (Л)=g;я+л-Зы (1+к 2)(1+д2к),

g;-т (Л) = g;_m - 2Л • ¿кт (1 + к2),

g-кт (Л) = g-кт - 2Л • ¿кт (1 + к2), £ к -т (Л) = £к-т + Л •¿т (1 + к 2)(1 + д ^ ),

¿кт - символ Кронекера,

М+

gкm = £ cosка{*Х) • cosта{]+1) + 1=1

м+

+ к • т •£ cosка(+2) • cosта(+2) + 1=1 ^ М1

+q

Ecos ka{ 1} ■ cos ma

j

(-1)+

j=1

+ k ■ m ■£ coska{'2) ■ cos ma{~2) ¿—i j j

j=1 Mj+

= £ cos kaj11 ■ si j=1

M+

sin ma()+l> +

+ k ■ m ■£ cos ka(+1) ■ sin ma(^2> +

j j

j=1 ^ MT

+q

k+m

Ecos ka1-. 1) ■ sin ma( 1) +,

j j

j=1

+ k ■ m ■£ coskaíJ-2) ■ sin?

v(-2)

J=1

g¡m = £ sin ka^+1) ■ cos ma^+1) + j=1

+ к • m •£ sin ка{.+2) • cos та(-+2) + ¿—i j j

j=1

(м+ -2 м- -2 ^

-к+м-)2 • Y\f(j>)\ +ц/л

+q

£ sin ка.'' • cos та, v + j=1

'+ к • m •£ sin ка{]2) • cos?

„(-2)

j=1

g'im = £ sin ка^1' • sin maj11 +,

+q

j=1

M+

к • m • £ sin ка<+2) • sin ma( ^ j j

j=1 ( M-

к +m

(+2)

Zsin ка( 1) • sin ma( 1) +

jj

j=1

(7)

Ag0-1 (Л) + Cg00№ + £ (0) CmgOm + £ (0) dmg0

m=-N m=-N

Ag0-lW + ^ Л) + £ (0) + £ (0) dmg0

m=-N m=-N

Agk.-1 + CÄ0 +£ (0) Cm§l(m Л) +£ (0) d„&m = ^ (к = ±1,..., ±N),

m=-N m=-N

Agh + CÄ0 +£ (0) Cmglm +£ (0) dmgm Л) = < (к = ±1,..., ±N)

м-

1 + к • m •£ sin ка(_2) • sin mа(f2) /—! j j

j=1

м+

м-

= £ фф +lnq•£ фф-q•£ у*,

j=1

м+ м-

e* = £ ф+ +£ ф

j=1 j=1

м2+

j=1

j=1

ek = £ • cos ка^1-1 + к £ уу • cos ка<^2) + j=1 j=1

^ м- м- N

+ q

£ фф • cos ка^ 1) - к • q •£ у *• cos ка^

-2)

v j=1 J- /

м1+ м+

e'k=£ ^^ • sin ка^1 + к £ уу • sin ка^2"1 + j=1 j=1 ^ м-

j=1

+ q

м-

£ ф- • sin ка^ 1) - к • q • £ у- * • sin ка{}

-2)

1 1=1 Матрица системы уравнений (7) является матрицей Грама, поэтому она положительно определённая и СЛАУ (7) имеет единственное решение.

4. Нахождение параметра регуляризации

Параметр регуляризации л> 0 найдём по методу невязки из уравнения

с(Л) = ¿0- (8)

Здесь с(Л) = с(А(Л), С0 (Л), с (Л), С-1 (Л),..., й-N (Л)), А(/), с (Л), с±1 (Л), й±1(Л),..., с± N (Л), й± N (Л) - решение СЛАУ (7) при фиксированном л> 0. Величину допустимой невязки ¿0 > 0 зададим с учётом требуемого уровня точности аппроксимации функций р± (£) и

У

(Z) на Г1 и Г2 :

j=1

j=1

(9)

Здесь ех, ех > 0, и е2 > 0 - заданные допустимые абсолютные погрешности аппроксимации функций р± у/± (£) на Г1 и Г2 соответственно.

Обычным образом [3] доказывается, что с(Л) -монотонно возрастающая функция при Л > 0 и с(Л) ^ при Л . Здесь

М1+ М[ м+ м-

с = £ (р()2+£ ()2 +£ ^*)2 +£ (у*)2.

1=1 1=1 1=1 1=1

Из сказанного выше следует, что уравнение (8) имеет (и притом единственное) решение Л = Л тогда и только тогда, когда выполняются неравенства С(0) ^¿0 <С„ . (I0)

В случае с(0) = ¿0 имеем Л = 0. Это соответствует случаю, когда коэффициенты с0, с±1, й+1, ...,с±к, по существу получены методом наименьших квадратов. При этом, как показывает практика, для больших номеров п коэффициенты с±п, й±п приобретают практически случайные значения даже при весьма малых погрешностях в задании функций р± (£), ц/± (О, определяющих краевые условия. Это

может привести к значительной неустойчивости решения - пульсациям И(£) на границе области Б и вблизи неё за счёт больших амплитуд высокочастотных гармоник функции Им .

Рассмотрим теперь случаи, когда условия (10) не выполнены. При с(0) > ¿0 уравнение (8) при данном

фиксированном номере N решения не имеет. Условие < ¿0 означает, что тождественно нулевая функция И(£) = 0 обеспечивает необходимый уровень допустимой погрешности аппроксимации граничных условий. Этот случай следует признать не представляющим практического интереса (заданная информация является некачественной).

Итак, в случаях, когда выполняется условие с(0) > ¿0, следует увеличивать N и повторять все

вычисления до тех пор, пока не станет с(0) < ¿0. Решение И(£) будет тем лучше в смысле требования о ^ т£, чем больше номер N при условии с < ¿0, так как при увеличении N величина

inf Q(c0

d_ N) на последовательности

^0' •••' - N '

расширяющихся множеств допустимых значений систем коэффициентов не увеличивается. Достигнув при каком-либо значении N требуемой точности приближённого равенства с = ¿0, увеличивать N далее нецелесообразно, так как увеличение затрат на возрастающий при этом объём вычислений не компенсируется повышением точности решения.

Решение уравнения (8) при выполнении неравенств (10) можно получить, применяя какой-либо численный итеративный метод, например, известный метод хорд.

5. Компьютерная реализация метода

По описанному алгоритму разработаны две программы: RingDirNeum-l.pas и RingDirNeum-2.pas, реализующие решение смешанной краевой задачи. Первая предназначена для построения конформного отображения f : D ^ K и служит вспомогательной

для второй, которая решает редуцированную краевую задачу для кругового кольца. Обе программы созданы в среде Turbo-Pascal.

Отметим их некоторые особенности.

1. В основу численного построения конформного отображения f : D ^ K положен итерационный метод альтернации [4-6].

2. Область D с границей Г моделируется в виде полигона DM (M = M1+ + M" + M + + M2- ), задаваемого множеством его вершин Z±ч Z±1),z(±2) Z(±2),

расположенных на Г . В случае необходимости осуществляется (в автоматическом режиме) сгущение заданного множества вершин путём его пополнения новыми (фиктивными) (с углами при них, равными п ), расположенными равномерно на прямолинейных отрезках, соединяющих пары исходных соседних вершин. Количество фиктивных вершин задаётся оператором произвольно по запросу программы. Производится (автоматически) пересчёт в каждую фиктивную вершину граничных значений h или dh/дп из двух ближайших смежных исходных - по закону линейной интерполяции. Цель описанного пополнения - повышение точности методов минимизации, используемых в процедуре альтернации при построении конформных отображений внутренности внешнего контура двусвязной области на внутренность круга и внешности внутреннего контура - на внешность круга.

3. На заключительном этапе построения отображения f : D ^ K, когда в качестве образа области D уже получена область, близкая к кольцу K, и последующие шаги итерации не приводят к улучшению отображения, используется метод малых граничных вариаций М.А. Лаврентьева [5, с. 25; 7].

4. На каждом шаге итерационных процедур (альтернации или граничных вариаций по Лаврентьеву) вычисляется модуль производной соответствующей отображающей аналитической функции во всех точках zj±2) g г± . Это необходимо для вычисления величин f '(Zf2))|, фигурирующих в формулах (3) и (9).

5. С помощью программ можно решать как смешанную краевую задачу, так и собственно задачи Дирихле и Неймана. В случае постановки задачи Дирихле вся граница Г области D описывается как Г+ и Г-(с заданием во всех граничных точках значений функ-

ции И), а в случае постановки задачи Неймана - как Г+ иГ- (с заданием значений дИ/ дп).

6. Опробование метода

Опробование метода проводилось обсчётом моделей двусвязных областей Р с различным набором граничных условий по двум вариантам: А) с точными граничными значениями известной гармонической функции и её нормальной производной; Б) с граничными значениями, осложнёнными случайной помехой с известными параметрами. Эти значения задавались в качестве исходных данных для программы Р1щР1т-Ывыш-1. Редуцированные граничные значения обрабатывались по программе Р1^Р1гЫвыш-2. Рассчитанные по формуле (4) в заданных внутренних точках области Р значения гармонической функции сопоставлялись с известными точными значениями, что давало возможность судить о точности достигнутых приближённых решений.

Приведём один пример опробования.

Модель «Эллипсы». Область Р ограничена эллипсами Г+ = {(х, у): х2 /25+ у2/9=1},

Г- = {(х,у): х2 /4+ у2 = 1} (рисунок). В качестве гармонической функции И взята функция И(х,у) = 0,5 + х - 0,2(х2 - у2) + 0,051п(х2 + у2) + 2,7(х + 3) 4,1(х - 9) _

(х + 3)2 + (у - 7)2 (х - 9)2 + (у + 2)2 - 0,4 х х2 + (у - 0,25)2 '

На границе Г области Р задано 180 точек, по 90 на каждой из граничных компонент Г+ и Г- . В половине из них на дугах (Л+ , В+), (Л+, В2+), (Л; , В; ), (Л,, В'2) заданы значения И, в остальных точках - значения дИ/дп. В качестве конформного образа области Р получено кольцо К = ^: 0,406 < Щ < 1}. При этом совершено 4 шага альтернации и 3 шага вариации граничных компонент по методу Лаврентьева. Максимальное

отклонение образа Г+ от единичной окружности составило 2,1-10-7, образа Г+ от внутренней окружности кольца - 2,8-Ю-4.

Расчёт по варианту А дал следующие результаты: максимальная абсолютная погрешность значений к(х, у), достигаемая в точке (-4,0;0,0), равна 0,07, что составляет 0,68 % от наибольшего значения вариации значений к (равной 11,4) в поле приведённой ниже таблицы. Среднеквадратичная погрешность значений к по всем 180 граничным точкам составляет 0,11, по всем внутренним точкам области Б в поле таблицы -0,02 (0,17 %).

Из таблицы видно, что максимальная абсолютная погрешность приближённых значений к(х,у), достигаемая в точке (-4,5;0,0), равна 0,25 (2,33 %). Достигнутое значение среднеквадратичной погрешности по всем внутренним точкам области Б в поле таблицы - 0,10 (0,93 %). В точках, удалённых от границы области на расстояние 1,0 и более, абсолютная погрешность составляет не более 0,17, причём погрешность уменьшается по мере удаления точки от границы внутрь области.

Анализ численных экспериментов позволяют сделать следующие выводы:

1. Точность метода достаточно высокая, особенно для случаев, когда граница области гладкая или (в случае области с негладкой границей) когда значения дк/дп задаются только на гладкой части границы.

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

В результате расчёта по варианту Б со случайной помехой, равномерно распределённой на промежутке (-0,10;0,10), при N =17, ^ = ^2 =0,2, получено

¿0 = 3,74 , Я = 0,1323. Приближённые значения к(х,у) во внутренних точках области Б в сопоставлении с точными значениями приведены в таблице.

Литература

1. Александров И.А. Конформные отображения односвязных и многосвязных областей. Томск, 1976.

2. Канторович Л. В., Крылов В. И. Приближённые методы высшего анализа. М., Л., 1962. С. 591.

3. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1979. С. 72.

4. Фильчаков П.Ф. Численные и графические методы математики: Справочник. Киев, 1970. С. 503-507.

5. Соболев В.В., Соболева (Ищенко) Н.В. // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2001. № 3. С. 23-26.

6. Соболев В.В., Соболева Н.В. Программа численного построения конформного отображения

Результаты решения смешанной краевой задачи Модель «Эллипсы» Значения к( х, у) (прибл. /точн.)

д -4,5 -4,0 -3,0 -2,5 -2,0 -1,5 -1,0 0,0 1,0 1,5 2,0 2,5 3,0 4,0 4,5

2,5 -0,90 -1,04 -0,06 -0,16 0,66 0,60 1,73 1,75 2,42 2,49 2,65 2,72 2,79 2,86

2,0 -3,38 -3,58 -2,32 -2,49 -1,38 -1,50 -0,55 -0,63 0,16 0,11 1,21 1,23 1,89 1,93 2,13 2,17 2,28 2,32 2,33 2,37 2,27 2,32

1,0 -6,39 -6,63 -4,00 -4,18 -2,95 -3,10 -2,01 -2,11 -1,16 -1,23 1,35 1,34 1,54 1,53 1,62 1,60 1,59 1,57 1,23 1,16

0,0 -7,93 -8,18 -6,60 -6,82 -4,22 -4,39 -3,17 -3,31 1,34 1,29 1,31 1,25 0,80 0,82 0,53 0,44

-1,0 -6,43 -6,62 -4,08 -4,21 -3,05 -3,15 -2,11 -2,20 -1,27 -1,36 1,29 1,26 1,47 1,41 1,53 1,46 1,48 1,39 1,03 0,93

-2,0 -3,62 -3,62 -2,55 -2,59 -1,63 -1,65 -0,10 -0,10 1,00 1,00 1,75 1,71 1,99 1,93 2,13 2,06 2,18 2,08 2,10 2,00

-2,5 -1,39 -1,21 -1,23 -1,21 -0,41 -0,38 0,30 0,34 1,44 1,47 2,24 2,19 2,48 2,41 2,62 2,53

ограниченной двусвязной области на круговое кольцо и обратного отображения. Ростов н/Д, 2001. Зарегистр. ГО ФАП ВНТИЦ, № 50200100349.

7. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М., 1965. С.377-378.

Ростовская государственная академия сельскохозяйственного машиностроения_5 мая 2006 г.