Об одном варианте метода прямых для уравнения Пуассона Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 12, № 2, 2007

ОБ ОДНОМ ВАРИАНТЕ МЕТОДА ПРЯМЫХ ДЛЯ УРАВНЕНИЯ ПУАССОНА*

А. М. Блохин Институт математики СО РАН им. С. Л. Соболева,

Новосибирск, Россия e-mail: Blokhin@math.nsc.ru

А. С. Ибрагимова, Н. Ю. Красников Новосибирский государственный университет, Россия e-mail: Ibragimova@ngs.ru

A theoretical justification of one algorithm for finding an approximate solution of the mixed boundary value problem for the Poisson equation is given.

Введение

В последние годы в вычислительной практике довольно широкое распространение получил метод прямых. Суть этого метода заключается в том, что в исходной дифференциальной задаче производится дискретизация только по части независимых переменных, т. е. исходное уравнение в частных производных аппроксимируется системой обыкновенных дифференциальных уравнений. Мы будем называть полученную таким образом вычислительную модель дифференциально-разностной моделью.

Однако возможны и другие приемы получения вычислительных моделей в методе прямых. Например, вместе с дискретизацией, скажем, по одной из независимых переменных можно использовать аппроксимацию производной по другой переменной с применением интерполяционных многочленов. Именно такой подход обсуждается в настоящей работе на примере уравнения Пуассона.

Поскольку при использовании метода прямых для конкретного случая мы сводим исходную проблему к краевой задаче для системы обыкновенных дифференциальных уравнений, возникает вопрос о нахождении ее приближенного решения. Хотя в настоящее время существует достаточно много алгоритмов численного решения краевых задач для систем обыкновенных дифференциальных уравнений, мы предлагаем в данной работе еще один способ нахождения приближенных решений таких краевых задач, использующий методы сплайн-функций.

Работа носит теоретический характер. Результаты численных экспериментов будут предметом для обсуждения в следующих работах.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 04-01-00900 и № 06-08-00384) и междисциплинарного интеграционного проекта фундаментальных исследований СО РАН-2006 № 46.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

1. Предварительные сведения

В настоящее время для описания физических процессов в полупроводниковых устройствах часто применяют так называемые гидродинамические модели (одна из моделей описана в работе [1]). Эта математическая модель содержит также уравнение Пуассона для электрического потенциала р (см., например, [1, 2]). На рис. 1 приведены граничные условия для потенциала р в случае полупроводникового транзистора с электронной проводимостью (подробно такой транзистор описан в работах [1, 2]). Итак, рассмотрим в области П уравнение Пуассона для потенциала р = р(х,у):

p = <fxx + ^yy = f (x, y)

(1.1)

с граничными условиями

py = 0 при y = 0, 0 < x < 6a; px = 0 при x = 0, 6a, 0 < y < 2a; p = 0 при y = 2a, 0 < x < a;

G = const при y = 2a, 2a < x < 4a; B = const при y = 2a, 5a < x < 6a; : 0 при y = 2a, a < x < 2a и 4a < x < 5a.

p p py

(1.2)

Здесь х,у — безразмерные независимые переменные; f (х,у) — известная правая часть; а > 0 — некоторая переменная. Не нарушая общности, положим далее, что а = п/6.

Далее, вместо (1.1) мы будем рассматривать параболическую регуляризацию этого уравнения:

pt = Ax,yp - f (x,y^ p = pC^^y^

(1.1')

где £ — независимая переменная, играющая роль времени. Краевые условия для уравнения (1.1 ) оставим теми же самыми (см. (1.2)). Добавляя к (1.1 ) и (1.2) начальное условие

p|t=o = po(x,y^

(1.3)

решение исходной задачи (1.1) и (1.2) будем искать, используя метод установления (по поводу этого метода см. работу [3]), т.е. вычисляя предел решения смешанной задачи (1.1'), (1.2) и (1.3) при £ ^ то.

Рис. 1. Область 0 и граница д0.

Как известно (см. [3]), для обоснования метода установления надо доказать асимптотическую устойчивость (по Ляпунову) стационарного решения ^(ж, у) смешанной задачи (1.1 ), (1.2) и (1.3). С этой целью введем функцию

и = и (Ь, ж, у) = ж, у) — -0(ж, у)

и сформулируем смешанную задачу для нахождения и:

и — Дх,уи = 0, Ь> 0, (ж, у) € П; (1.4)

и = 0 на Го, (1, Уи) = 0 на Гг, Ь > 0; (1.5)

и|4=о = ио(ж,у) = <£о(ж,у) — ^(ж,у), (ж,у) € П. (1.6)

Здесь

Г п п п 2 5

Го ^ (ж,у): у = з, 0 < ж < 6' з < ж < зп' 6п - ж - п

— часть границы дП (рис. 1); Гг = дП \ Го; 1 — единичный вектор внешней нормали;

Уи = (иж,и).

Умножим уравнение (1.4) слева на 2 и. После несложных выкладок мы получаем

(и2)4 — 2а1у(иуи) + 2|уи |2 = 0, |уи |2 = (уи, У и ) = их + и2.

Проинтегрируем полученное выражение по области П. Применяя формулу Остроградского — Гаусса, придем к следующему соотношению:

(1|и(Ь)||!2(п))4 + 2||Уи(Ь)||!2(п) — 2 У и(1, уи)^ = 0, (1.7)

да

где

||и(¿)|Ц2(п) = / и2(Ь,ж,у)^Ыу; п

||Уи(¿)|||2(П) = I |Уи(Ь,ж,у)|2^ж^у.

п

Используя известное неравенство Фридрихса [2]

||и(Ь)|^2(п) < С ||Уи(Ь)|ь2(п),

где С1 = С1 (П, Го) > 0 — постоянная, и учитывая краевые условия (1.5), из (1.7) получаем

2

(||и(¿)|||2(п)); + с||и(*)|Ц(п) < 0,

т. е.

_ 2 4

||и(^Н^п) < е ^ ||ио||2(п), Ь> 0. (1.8)

Здесь

||ио ^2(п) = ио2(ж,у)^ж^у.

п

Неравенство (1.8) означает, что решение -0(х,у) (решение исходной краевой задачи (1.1) и (1.2)) асимптотически устойчиво (по Ляпунову).

Для нахождения приближенного решения смешанной задачи (1.1 ), (1.2), (1.3) применим метод прямых. Производную рхх в уравнении (1.1 ) заменим не разностным отношением (как это обычно принято), а используем для аппроксимации ее интерполяционный многочлен. В соответствии с рекомендациями из [3] в качестве узлов интерполяции выбираем нули многочлена Чебышева. Пусть известны значения четной функции р(£,х, у) в узлах (рис. 2)

23-1 п, 3 = 17Ж

3 2Ж р(£,х3 ,у) = Рз (£,у).

Формула интерполирования четной функции р(£,х, у) будет иметь следующий вид [3]:

(1.9)

1

N

Р (х,Р) = (-1)1-1

,= 1

вт х,

сое х — сов х,

сов(Жх)р,(£, у), 0 < х < п.

(1.10)

Используя правило Лопиталя, из (1.10) получаем

Р(хз,р) = рз(£,у), з = 1,^.

Дифференцируя соотношение (1.10) по х нужное количество раз, придем к следующим выражениям:

дР(х, р) дх

1

ЛГ

N

^](-1)г 1 втх,р,(£,уК -

,=1

N вт(Жх) + сов(Жх) вт х

д2Р(х, р) _ 1

N

дх2

^](-1)г 1 втх,р,(£,уК -

сов х — сов х, (сов х — сов х,)2 N2 сов(Жх) N вт(Жх) вт х

(1.11)

,=1

2

сов х — сов х, (сов х — сов х,)2

+ сов(Жх) сов х +2 сов(Жх) вт2 х

(сов х — сов х,)2 (сов х — сов х,)3

(1.12)

Замечание 1.1. Из (1.11) следует, что

дР(х, р)

дх

0,

х=0, п

т.е. интерполяционная формула (1.10) автоматически удовлетворяет краевым условиям при х = 0 и х = п (см. граничные условия (1.2)).

Рис. 2. Узлы интерполяции.

Учитывая соотношения

cos(Nxj) = 0, sin(Nxj) = (—1)j+1, j = 1,N и используя правило Лопиталя, из (1.11) и (1.12) последовательно получаем

dp(x,

dx

N

£ (-D'+j-1 ■ÓnXw(lyl - 2ctgx,w(t,y),

Xj '=1, '=j

cos Xj — cos X' 2

d 2P (x,p)

dx2

N

= 2 £ (—i)

X=Xj '=1, '=j

N 2 - 1 1

'+j

-1 sin X' sin Xj^'(¿,y) (cos Xj — cos X')2

Замечание 1.2. Выражение преобразуем так:

3 — 2ctg2Xj ^j(t,y)-

2 sin X' sin Xj (cos Xj — cos X')2

(1.13)

(1.14)

2 sin X' sin Xj cos(2e'j) — cos(2a'j) 1f 1

(cos Xj — cos X')2 4 sin2 a'j sin2 Д? 2 [ sin2 ^'j sin2 a.

'j

где

i - j

Aj = о л т n, a.

i+j — 1

'j 2N ' 'j 2N

-n.

(1.15)

В уравнении (1.1 ) положим ж = Жj и агрегат заменим выражением (1.14). В

итоге будем иметь

= ^^ ^ ^ (Ь,у) + Е (—1) + /j (у), 3 (1.16)

'=1, '=j

Здесь a? = — + ^ — Й—i—; fj(y) = f (x?,y).

N2 111

--1-----

3 6 2 sin Xj

В (1.16) проведем также дискретизацию по переменной t. С этой целью введем следующие обозначения: ^j(y) = <£j(nA,y) = (y), n = 0,1,..., </3j(y) = (у), A _ шаг разностной сетки по переменной t.

. d^j (t,y) </?j (y) — (y) . . Аппроксимируя производную -—- выражением ---, перепишем (1.16)

так:

Здесь

dt

N

A

(y) = Aj&(y) + E (—1)'+jb'j&(y) + /(y), j = 1,n.

'=1, '=j

a=Aj+A, / (y)=/j (y)—

(1.17)

Систему (1.17) перепишем в векторном виде:

ф' = в Ф + т,

(1.17')

1

b

где

Ф

( р1(у)\ (у)/

т

/Л(у)\

\] N (у)/

В = ш + В; ш = diag(Дl,..., ДN) — диагональная матрица;

В

0

Ь 12

Ь 21 0

\Ь 1N Ь2N

ЬN Л ^^ 2

0

В*

Ь,з = (—1),+з Ь,з.

Замечание 1.3. Далее мы будем полагать, что т > 0 и матрица В = В* с доминирующей диагональю [4], т. е. В > 0 (за счет выбора шага А этого всегда можно добиться). К системе (1.17 ) надо добавить краевые условия (см. (1.2)):

Ф' (0) = 0,

Т2.

(1.18)

Здесь ^0>1 — диагональные матрицы порядка N (диагональные элементы равны либо нулю, либо 1), причем V) + = ^, ^ — единичная матрица порядка N; т — вектор размерности N (компоненты вектора т равны нулю, либо С, либо В).

2. Разрешимость краевой задачи (1.17') и (1.18)

Обсудим теперь очень важный для всей работы вопрос о разрешимости сформулированной выше краевой задачи (1.17) и (1.18). Поскольку матрица В = В* > 0, рассмотрим матрицу Л = В1/2. Матрицу Л можно определить так. Существует ортогональная матрица Т = Т(В) такая, что [4]

В_Т^^Т* т—1_Т*

где Д = diag(A1 ,...,ЛN) — диагональная матрица порядка N; Лз- = Лз-(В) > 0, з =

— собственные значения матрицы В. Тогда Л = В1/2 = ТД1/2Т*. Здесь diag(^/Al, ..., ).

Переформулируем краевую задачу (1.17) и (1.18) так. Сконструируем вектор

^(у) = (ф

где Ф = Ф — ЛФ. Тогда систему (1.17 ) можем переписать следующим образом:

V' (у) = МУ(у) + Р(у). (2.1)

Здесь М = ( л ^ ) — квадратная блочная матрица порядка 2N; ON — квадратная \ON —Л/

нулет.я п°рядка N; р = — векТОр _ частей ра3мерности ^; 0N —

нулевой вектор размерности N.

Краевые условия (1.18) примут такой вид:

(Л,/„ )У(0) =0,

(*1,*)у( §)= <2'2)

где (Л, /^), (/>1, ^о) — блочные матрицы размерности N х 2Ж; />1 = ^оЛ +

Решение краевой задачи (2.1) и (2.2) ищем в таком виде (см. по этому поводу также [5]):

y

' -s)

V/(y) = eyMV(0) + J e(y-s)MF(s)ds, 0

причем вектор V(0) определяется из системы линейных алгебраических уравнений

KV (0) = /• (2.3)

Матрицу K и вектор / выпишем ниже. Сначала же покажем, как вычисляется матричная экспонента eyM. С этой целью воспользуемся так называемым разложением Холес-ского [6]

M = Q block diag(Л, -Л)^-1, (2.4)

где

q i x) q-1 _ ( -х

q l.ow /J ' q vow in

блочные квадратные матрицы порядка 2N; х _ -1/2 ■ Л 1. Тогда (см. [5, 6])

feyA еуЛЛ-1 - Л-1е-уЛ\

eyM _ Q block diag(eyV-^)Q-1 _

блочная квадратная матрица порядка 2N, при этом еуЛ _ Tdiag(eyA,/^1 , ...,вУл,/^)T*

\On е-уЛ

ey

квадратная матрица порядка N.

Вычисляя вектор У и подставляя его в (2.2), выпишем матрицу K и вектор 0

Л In

K = |_ пЛ „eПЛЛ-1 - Л-1е-3Л Л i/ie зЛ i/1----+ v0e зЛ

квадратная блочная матрица порядка 2N; 0 = ( ) — вектор размерности 2N;

i J

>2 _ V2 - (/>1, vo)e^ J e"*MF(s)ds.

0

Система (2.3) однозначно разрешима, если

det K _0. (2.5)

3

Поскольку

е 3 ^л—1 _ Л—1е—3л п п

det К = det Л det <| ^--+ ^0е-3л — ¿^е3лЛ-

е 3 л Л-1 + Л-1е—3

det Л det — ^-^--± ^0е

,-П л

3=

=—3л _ е3л 3л , „3л = det Л det т/0--г^Л'

— — л — л — — л I — л

е 3л — е3л 4 _1 е 3л + е3' ч 2"

— 0 det Л det Т* det |т/0Тdiag(1 — е-3,..., 1 — е-3) + / 1 ± е-1плАГ 1 ± е-1^лА^Ч 1

+-Ц 1±кН Ь (е3 Д) =

1 ^ N /1 ±е-2п^\ * detЛdetТ* det (е3л) detТП^ уА~-)Х

х det щТdiag( л/А-,..., \fANN-2 )Т* ± щ

[ \ 1 ± е- 31 ± е-3)

условие (2.5) выполнено, если

det (щ -А ± ^) = 0. (2.5')

Здесь

/ _1 _ е-1п\^Л1 ._1 _ е-3

А = А* = Т^ л/А1-2ГТТ, ...^л/А^-^^ Т* > 0.

V 1 ± е-3^Л1 1 ± е-2^VЛN у

Рис. 3. Матрица ^0 - А + ^1.

На рис. 3 схематично изображена матрица •А + VI. Здесь заштрихованы массивы матрицы ^ •А + VI, равные соответствующим массивам матрицы А. Следовательно,

det •А + VI) = det а • det а2 = 0,

поскольку А > 0. Таким образом, решение краевой задачи (1.10') и (1.11) окончательно будет иметь следующий вид:

у

V (у) = еуМ •К-1 • 0 + [ е(у_*)м • / (з)^. (2.6)

3. Нахождение приближенного решения краевой задачи (1.17') и (1.18) с помощью техники сплайн-функций

Учитывая трудности, возникающие при практическом использовании формулы (2.6) для нахождения приближенного решения краевой задачи (1.17 '), (1.18), мы предложим следующую технологию построения ее приближенного решения. Будем искать приближенное решение в виде интерполяционного кубического сплайна класса С2 [7]:

й2

Б(у) = (1 - т)Фк + тФ^1 - -Ут(1 - т)[(2 - т)т + (1 + т)тк+1], (3.1)

6

у _ Ук П ''

где т = —й—, у € [Ук,Ук+1 ], Ук = кйу, к = 0,К - 1, Кйу = 3; Фк = Ф(у*.); т = Ф"(у*.).

Кубический сплайн (3.1) должен быть непрерывным вместе со своей первой производной всюду на отрезке [0,п/3] (вторая производная непрерывна по определению сплайна (3.1)). Вычисляя векторы

$(ук + 0), $(ук - 0)

и приравнивая их, получим

1 13 _

2Шк_1 + 2т*, + 2= (пФк - ПФк), к = 1, К - 1. (3.2)

Здесь п = Ф - 1, П=1 - Ф-1 — разностные операторы; ф, ф-1 — операторы сдвига: ф^Ф* = Фк±1, Ф+1 = ф. " '

Полагая в (1.17 ) у = у^ и подставляя Ф (у&) (см. уравнение (1.17 )) в (3.2), получим следующую разностную схему:

- й6УФк_1 - 21/* + уя}Фк + { ^ - й6Уя}Фк+1 =

й2 ^ ^ ^ _

= -У+ 4. + к = 1,К - 1, (3.3)

6

где . = .(ук). _

Для нахождения сеточной вектор-функции Ф&, к = 0,К, с помощью системы алгебраических уравнений (3.3) необходимо задать краевые условия при к = 0 и к = К.

Аппроксимируя в краевых условиях (1.18) производную разностным отношением, в итоге получаем

Ф0 = Ф1,

0 - } (3.4)

Фк = £Фк-1 ± Ь. /

Здесь С — диагональная матрица порядка N (диагональные элементы равны либо нулю, либо 1); Ь — вектор размерности N (компоненты вектора Ь равны либо нулю, либо С, либо В). Систему алгебраических уравнений (3.3) и (3.4) можно решить методом матричной прогонки [8]. Легко видеть, что из неравенства В > 0 следуют условия хорошей обусловленности разностной краевой задачи (3.3) и (3.4).

Список литературы

[1] Romano V. 2D simulation of a silicon MESFET with a non-parabolic hydrodynamical model based on the maximum entropy principle //J. Comput. Phys. 2002. Vol. 176. P. 70-92.

[2] Blokhin A.M., Bushmanoy R.S., Rudometoya A.S., Romano V. Linear asymptotic stability of the equilibrium state for the 2-D MEP hydrodynamical model of charge transport in semiconductors // Nonlinear Analysis. 2006. Vol. 65. P. 1018-1038.

[3] Бабенко К.И. Основы численного анализа. М.; Ижевск: НИЦ "Регулярная и хаотическая динамика", 2002.

[4] Ланкастер П. Теория матриц. М.: Наука, 1980.

[5] Годунов С.К. Матричная экспонента, матрица Грина и условие Лопатинского: Учеб. пособие. Новосибирск: Изд-во НГУ, 1983.

[6] Блохин А.М. Элементы теории гиперболических систем и уравнений: Учеб. пособие. Новосибирск: Изд-во НГУ, 1995.

[7] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980.

[8] Самарский А.А., Попов Ю.П. Разностные схемы газовой динамики. М.: Наука, 1975.

Поступила в редакцию 14 сентября 2006 г.