CC BY
48
МАТЕМАТИЧЕСКАЯ ФИЗИКА, МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
MSC 65N12
АНАЛИТИКО-ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ОДНОГО КЛАССА СИНГУЛЯРНО ВОЗМУЩЕННЫХ СИСТЕМ НЕЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
С.И. Безродных* **, В.И. Власов*
*Вычислительный центр им. А.А.Дородницына РАН, ул. Вавилова, 40, Москва, 119333, Россия.
**Государственный астрономический институт им. П.К. Штернберга МГУ, Университетский пр., 13, Москва, 119992, Россия, e-mail: sergeyib@pochta.ru, vlasov@ccas.ru
Аннотация. Для одного класса сингулярно возмущенных систем нелинейных дифференциальных уравнений предложен высокоэффективный аналитико-численный метод решения, основанный на сочетании операторного метода Ньютона и метода продолжения по параметру. а также на построении начального приближения в явном виде. Эффективность метода продемонстрирована на важном частном случае этой системы. возникающем при моделировании взаимодействия физических полей в полупроводниковом диоде. Для этого случая производная Фреше и функция Грина соответствующего дифференциального уравнения найдены аналитически. Численная реализация подтвердила высокую эффективность и сверхэкспоненциальную скорость сходимости предложенного метода.
Ключевые слова: сингулярно возмущенные системы, аналитико-численные методы, моделирование полей.
1. Введение
При моделировании ряда важных физических явлений (в том числе включающих процессы перноса, диффузии и дрейфа заряженных частиц различных типов [1]- [6]) возникает краевая задача для сингулярно возмущенной системы из M > 3 дифференциальных уравнений следующего вида:
[e2(x) Ф'(ж)]' + Ф(я) , P(x)) = signx , (1)
\ат (Ф) Р'тО) + Ф'О) Ът (Ф) Рт О) ]' = Пт (х, Р), т = 1, М - 1, (2)
рассматривать которую будем на отрезке [-1,1]. Искомыми величинами здесь являются скалярная функция Ф(х) и вектор-функция P(x) = |Pi(x),... , PM-i(x)}.
Работа выполнена при финансовой поддержке РФФИ (проект №13-01-00923), Программы ОМН РАН «Современные проблемы теоретической математики», проект «Оптимальные алгоритмы решения задач математической физики» и Программы №3 фундаментальных исследований ОМН РАН.
НАУЧНЫЕ ВЕДОМОСТИ Е|Я Серия: Математика. Физика. 2013. №19(162). Вып. 32 149 Для системы (1), (2) рассматриваются граничные условия Дирихле:
Фигурирующая в уравнении (1) величина е(х) представляет собой ступенчатую функцию, определяемую по формуле
где е± — малые положительные параметры, а присутствующая в том же уравнении величина £(Ф, Р) является гладкой функцией от Ф и Рт, для которой выполняются условия согласования
уравнения (1) с краевыми условиями (3). В уравнениях (2) коэффициенты am, bm являются гладкими функциями от Ф, причем am не обращаются в нуль, а правые части Rm (x, Р) этих уравнений достаточно гладко зависят от аргументов x и Pm.
Уравнения (1), (2) понимаются в смысле выполнения соответствующих интегральных тождеств, что эквивалентно их поточечному удовлетворению на (-1, 1) \ {0} и выполнению известных условий трансмиссии в точке x = 0. При этом предполагается, что функции Ф и Pm принадлежат классу B := C2(I-) П C2(I+) П C[— 1, 1], где I± определяются равенствами I- := [-1, 0], I+ := [0, 1].
В работе получены следующие результаты. Предложен аналитико-численный метод решения краевой задачи (1)-(3), основанный на сочетании операторного метода Ньютона [7]- [9] и метода продолжения по параметру [10], а также на построении начального приближения в явном виде, см. разд. 2. Высокая эффективность метода продемонстрирована на примере задачи, возникающей при моделировании физических полей в полупроводниковом диоде, см. разд. 3, 4.
Полученные результаты является развитием [11]- [18].
2.1. Сведение к системе интегро-дифференциального и (М —1) интегральных уравнений.
Обозначим через Ят(х) правые части уравнений (2), понимаемые как функции х,
и перейдем от исходной задачи, сформулированной относительно неизвестных {Ф, Р}, к задаче для неизвестных {Ф, И.}, где И. — вектор, образованный величинами (4):
Заметим, что дифференциальные операторы в левых частях уравнений (2) линейны относительно искомых Рт и имеют специальный вид, который позволяет выписать
Ф(±1) = ±V , Рт(±1) = Р,^ , т = 1, М — 1.
(3)
x € [-1, 0) , x € (0,1] ,
£(± V, Р(±1)) = ±1,
1. Метод решения краевой задачи
Rm (#) = Пт(х-, Р) , т = 1, М — 1 ,
(4)
(5)
явно функции Грина РШ(Ф; х, £), соответствующие однородным условиям Дирихле на концах отрезка [-1,1]:
Р (Ф; ±1,0 = 0
Эти функции Грина даются равенствами
Рт(Ф; х,£) := - е-1т(х) Ел1
X ; а-1
ХШ) ит
Е
Хш; аш
Е
X ; а-1
ХШ) 'Ш
(6)
Рш(Ф; х,е) = РШ(Ф; х,С), С е (-1,1), х е I±К), (7)
где 1-(е) := [-1, £], 1+(е) := [С, 1]) а функции РШ определяются выражениями
(8)
РШ(Ф; х, С) := - е
Здесь использованы обозначения
гв
-1
X ; а-1
Хш) иш
Е
Хш; аш
Е
X ; а-1
хт аш
(9)
Е^ [и; V] := / г>(£) ехр [«(£)] ^ , Е [и; V] := Е-1 [и; V]
-1
1
аш
— . (ю)
а ХШ(х) определяется через коэффициенты уравнения (2) по следующей формуле:
Ф(ж)
1„,М:= / М§-<Й,
V
аш(^)
т = 1, М — 1.
(11)
Используя функции Грина РШ, определяемые из (7)-(9), и вводя обозначение для свертки
<#(х, е) * /(е) > := ^1 #(х, е) /(е)
получаем из уравнений (2) с учетом граничных условий (3) и обозначений (4), (5) следующие выражения для функций РШ через функцию Ф и вектор-функцию И.:
Рт(х) = Тт [Ф, II] , т = 1, М - 1;
Тт [*, П] := (Р„,(*; *. О* Д,.,«)) + е~1тМ { П. + (р^ | .
(12)
(13)
Первые слагаемые в правых частях равенств (13) представляют собой решения уравнений (2) с однородными граничными условиями, а вторые — решения соответствующих однородных уравнений с неоднородными граничными условиями (3).
Заметим, что правые части формул (8), (9), (13) зависят от Ф через входящие в них коэффициенты аШ уравнения (2), а также величины Хш, определяемые из (11).
Подставляя Рт = Vm [Ф, R], т = 1 ,М— 1, из формул (12), (13), в уравнение (1), приходим к нелинейному интегро-дифференциальному уравнению относительно Ф, которое вместе с граничными условиями из (1) составляет краевую задачу
D (R) [Ф] = 0 , Ф(-1) = -V, Ф(1) = V; (14)
здесь интегро-дифференциальный оператор D (R) [Ф], действующий на функцию Ф(х), определяется по формуле
D (R) [Ф] := [е2(ж)Ф'(ж)]/ + L(R)^] — sign x . (15)
где введены обозначения
L(R)[Ф] := L(Ф(х), P [Ф, R]) , Р[Ф, R] := {Pi^,Ri],..., Pm-i^,Rm-i]} . (16)
Отметим, что для оператора L и уравнения (14) функция R(x) играет роль параметра.
Подставляя Pm = Pm [Ф, R] из формул (12), (13) в (4), получаем (M — 1) интегральное уравнение для функций Ri(x),... RM-i(x):
К-m (Ф)[Р] := Rm - ftm(*,P [Ф, R]) = 0, 1П (17)
интегральные операторы Km (Ф)[R], определяемые равенствами (17), действуют на функцию R(x), а функция Ф(х) играет роль параметра.
Таким образом, задача (1)-(3) для функций {Ф, Р}у сведена к краевой задаче для интегро-дифференциального уравнения (14) и интегральных уравнений (17) относительно функций {Ф, R}y. Индекс V в приведенных обозначениях подчеркивает зависимость решений задачи (1)-(3) или задачи (14), (17) от параметра V из граничных условий для Ф(х).
2.2. Итерационный алгоритм. Для решения задачи (14), (17) предлагается метод, основанный на операторном методе Ньютона [7]- [9] и эффективном построении начального приближения для него при помощи способа, описанного ниже. Указанный итерационный процесс записывается в следующем виде:
Фга+1 = Фп — Р(Фга, Rn) [D(Rn)^ra]l , (18)
Rnml = Кг - Кт(Ф'\ R") [/Ст(Фга) [1C]], т = 1,М-1 , (19)
где верхний индекс, означающий номер итерации, принимает значения n = 0, 1, ...
В формулах (18), (19) линейные операторы D и Km действуют на аргументы в квадрат-
ных скобках и параметрически зависят от аргументов в круглых скобках. Они являются обратными операторами к производной Фреше соответственно для операторов D и Km.
Таким образом, приближенное решение задачи (14), (17) строится по формулам (18), (19), а ее точное решение ищется в виде предела {Ф, R}V = lim {Фп, Rn}V. Подставляя
предельные функции Ф(х) и R(x) в формулы (12), (13), вычисляем Рт(х), т = 1, М — 1, чем и завершаем решение исходной задачи (1)-(3).
2.3. Построение начального приближения. Построение начального приближения {Ф0, R0} V для итерационного процесса (18), (19) представляет собой самостоятельную трудную задачу. Предлагаемый аналитический способ построения такого приближения (см. [15], [17]) заключается в том, чтобы в качестве {Ф0, R0}v выбрать функции Ф0 и R0 = R(x, Р0) = {R1(x, Р0),... , RM(x, Р0)}, где Ф0 и Р0 — решение модифицированной системы (1), (2), в которой первое уравнение оставлено без изменения, т.е.
Ге2(х) -^-Ф°(х)1 + с(ч>°(х), Р°(х)^ = signer, (20)
dx L dx J V /
а остальные (M — 1) уравнений упрощены, например, исходя из физического смысла системы (1)-(3). Отметим, что такой подход позволяет построить {Ф0, R0}v как правило для ограниченного диапазона параметра V, фигурирующего в граничных условиях (3). Для других значений V начальное приближение строится при помощи метода продолжения по параметру [10], заключающегося в следующем.
Пусть для V = V* начальное приближение ^0,R0}v* известно, а требуется найти решение задачи (1)-(3) при V = V*. Согласно сказанному решение ^,R}V* строится при помощи процедуры (18), (19).
Чтобы получить начальное приближение при V = V* разобьем отрезок [V*, V*] точками Vk (где k = 0, K) так, что
V* =: V0 < Vi < ...Vk < ... < Vk := V* , (21)
и обозначим через {Ф, R}Vk решение задачи (14), (17), соответствующее V = Vk и состоящее из функций Ф^), R(Vk). Здесь и далее аргумент Vk, указанный в круглых скобках, означает, что функция соотсетствует V = Vk в граничном условии (1) или (14). Если решение известно для V = Vk-i, то в качестве начального приближения {Ф0, R°}Vfc для {Ф, R}Vk выбираем функции
Ф0^) := Ф^-i) + (Vk — Vk-i) x , R0(Vk) := R(Vk-i); k =1, 2,...,K , (22)
и при помощи итерационного процесса (18), (19) получим в пределе
Ф^) = lim Фп^), R(Vk) = lim Rn(Vk). (23)
Поскольку решение {Ф, R}v0 известно, то процедура (23) с начальным приближением (22) дает продолжение решения от меньших k к большим. При k = K с помощью (22), (23) получаем начальное приближение {Ф0, R0}v* и само решение {Ф, R}v* задачи (14), (17). Подставляя предельные функции Ф^*), R(V*) в формулы (12), (13), завершаем построение решения {Ф, Р^* задачи (1)-(3).
3. Приложение к расчету полей в полупроводниковом диоде
3.1. Постановка задачи. При моделировании полупроводникового диода в одномерном диффузионно-дрейфовом приближении [19]- [23] возникает частный случай рассмотренной выше краевой задачи (1)-(3). При этом функция Ф является безразмерным электрическим потенциалом. Вектор Р состоит из двух функций, обозначаемых P и N,
Р(х) = {P(x), N(х)} ,
которые представляют собой соответственно безразмерные плотности дырок и электронов. Функция sign x, фигурирующая в правой части (1), описывает плотность «вмороженных» зарядов. Будем рассматривать случай, когда материал диода однородный и, следовательно, е не зависит от х; эта величина представляет собой малый параметр порядка 10-4. Для рассматриваемого случая, возникающего при моделировании полупроводникового диода, функции L, am, bm, Rm, фигурирующие в системе (1), (2), имеют следующий вид:
L(P, N) = P — N , ац(Ф) = а2(Ф) = МФ) = 1, ЫФ) = —1, (24)
Ri(x, R) = R, R2(x, R) = R,
где и суть малые параметры порядка 10-4; R(x) — рекомбинационная функция, взятая в виде, предложенном Шокли, Ридом и Холлом [1],
г PN-N2
R(x) = П Р, N := ------------, (25)
V ; L J N + Nt + T{P + Nt)' y 1
где параметр т лежит в диапазоне [0.1,10], а параметр N имеет порядок 10-9. Из (24) следует, что функции Ii и I2, определяемые в (11) принимают в данном случае вид
Ii(x) = Ф(х) — V, I2 (х) = Ф(х) + V. (26)
Отметим еще, что величина V £ [—100,100], фигурирующая в граничных условиях (3), является половиной безразмерной разности потенциалов, подаваемой на диод; граничные условия для концентраций дырок и электронов имеют вид: P(—1) = N(1) = 0.
P(1) = N(—1) = 1. Таким образом, для рассматриваемого случая задача (1)-(3) пре-
вращается в следующую (см. [14]):
е2 Ф// + P — N = sign x , (27)
P// + Ф^ + Ф^ = £P R , N// — Ф^ — Ф^ = R; (28)
Ф(±1) = ±V ; P(—1) = 0, P(1) = 1; N(—1) = 1, N(1) = 0 . (29)
Решение задачи (27)-(29) ищется в классе £> := C2[—1, 0] П C2[0,1] П Ci[— 1, 1].
Задача (27)-(29) является сингулярно возмущенной, поскольку множителем при старшей производной в уравнении (27) служит малый параметр е2 .
3.2. Сведение к системе двух уравнений. Поскольку краевая задача (27)-(29) является частным случаем задачи (1)-(3), то в соответствии со сказанным в разд. 2 она сводится к задаче, аналогичной (14), (17), где вместо (14) фигурируют следующие соотношения:
D (R) [Ф] := е2 Ф//(х) — sign x + L (R)^], Ф(—1) = — V , Ф(1) = V, (30)
а вместо системы (17), состоящей из (M — 1)-го интегрального уравнения, возникает одно интегральное уравнение
K (Ф)[^ := R — Rp [Ф, R], N [Ф, R]j =0; (31)
где R выражается по формуле (25).
Интегральный оператор L (R) [Ф] в (30) определяется выражением
С (Я) [Ф] := ev-* - е1'+* + 6р { Р(Ф; .г, £)*/?«)> - In < Г»(Ф; *,
(32)
где использованы обозначения
Г в
Ef [u] := / exp [u(t)] dt, E [u] := E-i [u] . (33)
J a
Функции Грина в свертках (32) даются формулами
Р(Ф; х, С) = Р±(Ф; х, С) , Н(Ф; х, С) = Н±(Ф; х, С) , x £ 1±(С) , в которых для Р± и N± справедливы формулы, вытекающие из (8), (9):
Р-(Ф;.г,С) := , Р+(Ф;.г,С) := [Ф] -^21;
N±(Ф; х, С) = Р±(—Ф; х, С).
Для функций P(х) и N(х) имеем вместо (12), (13) следующие представления через Ф и R:
Р(х) = V [*, В] := 6р (Р(Ф; х, О . Д(0> + ev-*w -^1, (34)
ЛГ(*) = М\р,Щ fH <N(*:*.0 . ЩО) + ev+*w ^Г.^1 . (35)
Для решения задачи (26)-(28) в рассматриваемом случае, как и в общем, применяется изложенный в разд. 2 метод, в сочетании с эффективным построением начального приближения. Итерационный процесс (18), (19) для данного случая имеет вид:
Фп+1 = фп — 0(фп, Rn) [D(Rn)^n]l , (36)
дп+1 = дп _ ^(фп, д") К(Ф")[Дп]
(37)
Операторы О и К, фигурирующие в (36), (37), построены в аналитическом виде, см. разд. 4, а начальное приближение {Ф°, Я0} для итерационного процесса взято из [15], [17]. Заметим также, что при реализации процесса (36), (37) удалось избежать численного дифференцирования.
4. Численная реализация для модели полупроводникового диода
4.1. Производные Фреше операторов Р и К. Для того чтобы построить операторы Ь и К, фигурирующие в процессе Ньютона (36), (37), приведем вначале вид обратных к ним, т.е. производных Фреше операторов Р и К. Производная Фреше Р'(м,/) [г>] оператора Р (/) [и] из (30), где £(/)[и] дается равенством (32), имеет вид
V' (и, f) [^] = є2 ^"(ж) _ Т (и, f; ж) V (ж) + 2(м, f) [V]
(38)
а производная Фреше К' (и, f) [V] оператора К (и) ^] из формулы (31), где К определяется из (25), имеет вид
К' («. f) (V] = V (ж) _ <5,, Л (и, f) [V];
(39)
выражения для Т, Q и Л не приводятся здесь в силу их громоздкости.
4.2. Обращение производной Фреше оператора V. Обозначим через Фп+1 разность функций Ф"+1 и Фп, т.е. фп+1 := ф"+1 _ фп. Построение оператора О, фигурирующего в процессе Ньютона (36), эквивалентно решению следующей краевой задачи относительно ф п+1: ^ ^
V'(Фn, Д") [Фп+1] = ^(Я") [Фп] ; Фп+1 (±1) = 0 . (40)
Функции фп+1 будем искать в виде ряда
ф"+1 = £ х"+
т=0
(т)
члены которого хП+1 являются решением следующих краевых задач
Н(Ф", Я") [х^] = -Р(Я") [Ф"]; хП°+1 (±1) = 0. Н(Ф", Я") [х"?!] = -е(Ф", Я") [х"+-11)]; х"+1 (±1) = 0, т =1, 2,
где оператор Н (и, /) [V] дается равенством
Н (и, f) [V] := є2 v''(ж) _ Т (и, f; ж) V (ж)
Решения х(+)1 задач (42), (43) представляем в виде свертки
Х'"+. = _(С(Ф",Я";) * V(Я") [Ф"])
(41)
(42)
(43)
(44)
х"+ = _(С(Ф", Д"; ж, 0 * Q(Ф", д") [х"+Т1) ]) . т =1. 2,..., (45)
правых частей уравнений (42), (43) с функцией Грина С(и, f; ж, £) оператора Н (и, f) [V] с однородными условиями Дирихле v(±1) = 0. Таким образом, фигурирующий в (19) оператор О, обратный к производной Фреше оператора V, находим в виде
О(Ф", дп) [р(Д")[Ф"]1 = _ф"+1, (46)
где функция ф"+1 дается формулой (41).
Поскольку параметр є мал, то для вычислительных целей достаточно использовать вместо функции С главный член С0 ее асимптотики при є ^ 0, вид которого, полученный методом ВКБ [24], следующий:
С- (ж, О = _є-1 [Т(ж) Т«)Г1/4 Л (Л) Л (^£‘) / * (.7) , (47)
г/3
<#(*. о = €£({,*). Я ■■=€-' /гЩМ, 3-.= 3±^ (48)
о а
в этих формулах для краткости С±(и^; ж, £) заменено на С±(ж, £), а Т (и, f; ж) — на Т (ж).
4.3. Обращение производной Фреше оператора К. Обозначим через Дп+1 разность функций Дп+1 и Д", т.е. Дп+1 := Дп+1 _ Д". Построение оператора К, фигурирующего в процессе Ньютона (37), эквивалентно решению следующего уравнения относительно Д"+1:
К' (Ф", Д") [Д^1] = _К (Ф") [Д"].
Функции Д"+1 будем искать в виде ряда
Дї"+1 = £ Лт п"+!, (49)
"р
т=0
в котором функции п"+1 даются следующими формулами:
ч"+1 = _К (Ф") [Д"]; ч"’+! = Л (Ф", Д") [,"+-‘)] , т = 1, 2,...
Таким образом, оператор К из (37), обратный к оператору производной Фреше оператора К, находим в виде
К(Ф", Д") ІК(Ф") [Д"]1 = _-ф"+1 , (50)
где функция Д?"+1 дается формулой (49).
4.4. Реализация итерационного процесса. На нулевом шаге итерационного алгоритма строится начальное приближение {Ф°, Я°}у способом, описанным в п. 2.3. Точнее говоря, при V < у < 0 взято начальное приближение, построенное [15] в аналитическом виде, а для произвольных V приближение строится методом продолжения по параметру V. На первом и следующих шагах алгоритма вычисляются функции Ф" и Я" при помощи процесса Ньютона (36), (37), где для О и К используются формулы (46) и (50). Решение задачи (14), (31) ищется в виде предела {Ф, Я}у = Иш {Ф", Я"}у, после чего
функции Р и N из (27)-(29) находятся подстановкой вычисленных функций Ф и Я в формулы (34) и (35). Таким образом, решение {Ф, Р, N}у задачи (27)-(29) построено. О сходимости итерационного процесса см. в [17].
4.5. Замечание об эффективности метода. Изложенный метод решения задачи (27)-(29) был численно реализован. Искомые функции были вычислены для широкого диапазона изменения входных параметров. Приведённые в работе численные данные соответствуют следующим их значениям:
є = 4.35 ■ 10-4 , ЛР = 1.731 ■ 10-4 , = 0.577 ■ 10-4 , т = 1, N = 10-9 .
Результаты, представленные на рис. 1-4, получены при V = 0, на рис. 5-8 — при V = 6, Численные эксперименты показали, что для указанных значений входных параметров можно положить ф = _1.
На рис. 1, 5 даны графики функции Ф(ж) на всём отрезке [_1, 1], на рис. 2, 6 — графики той же функции на отрезке [_0.004, 0.004], т.е. вблизи начала координат, что позволяет более детально представить структуру внутреннего слоя; на рис. 3, 7 даны графики функции Р (ж) на отрезке [_1, 1]; на рис. 4, 8 даны графики функции N (ж) на отрезке [_1, 1].
Метод показал высокую эффективность при всех рассматривавшихся значениях входных параметров. В частности, данные, представленные на рис. 1-8, получены с относительной погрешностью 10-15 при использовании 5 итераций алгоритма (36), (37).
Сверхэкспоненциальная скорость сходимости наблюдалась при всех рассматривавшихся значениях входных параметров и проиллюстрирована на рис. 9, 10. На рис. 9 изображен график приближенного решения Ф4(ж) задачи (27)-(29) при V = 1. На рис. 10 изображены графики функций Ф”(ж) = Ф” — Ф”-1, п =1,4. Графики демонстрируют следующее соотношение |Ф" _ Ф"-1| = С ехр ( _ 1.73 п2), означающее указанную скорость сходимости.
Рис. 1. График функции Ф(ж) на отрезке [-1,1] при V = 0.
Рис. 2. График функции Ф(ж) на отрезке [-0.004, 0.004] при V = 0.
Рис. 4. График функции N (ж) на отрезке [-1,1] при V = 0.
Рис. 5. График функции Ф(ж) на отрезке [-1,1] при V = 6.
Рис. 7. График функции Р(ж) на отрезке [-1,1] при V = 6.
Рис. 8. График функции N (ж) на отрезке [-1,1] при V = 6.
Ф1 (ж)
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
-0.1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
X
(а)
Ф2(ж)
(б)
(в) (г)
Рис. 10. Иллюстрация сверхэкспоненциальной скорости сходимости метода при V = 1.
Литература
1. Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред / М.: Гостехиздат, 1957. 2-е изд.: М.: Наука, 1982.
2. Куликовский А.Г., Любимов Г.А. Магнитная гидродинамика / М.: Физматгиз, 1962. 2-е изд.: М.: Логос, 2005.
3. Половин Р.В., Демуцкий В.П. Основы магнитной гидродинамики / М.: Энергоатомиздат. 1987.
4. Морозов А.И. Введение в плазмодинамику / М.: Физматлит, 2006.
5. Gaevsky Kh., Groger K. and Zakharias K. Non-linear operator equations and operator differential equations / Mir: Moscow, 1978.
6. Agoshkov V.I., Dubovski P.B., Shutyaev V.P. Methods for solving mathematical physics problems / Cambridge International Science Publishing, 2006.
7. Канторович Л.В., Акилов Г.П. Функциональный анализ / М.: Наука, 1977.
8. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы / М.: Наука, 1987.
9. Колмогоров А.Н., Фомин С.В. Элементы теории функций и функционального анализа / М.: Наука, 1976.
10. Федоренко Р.П. Введение в вычислительную физику / М.: МФТИ, 1994.
11. Vlasov V., Bezrodnykh S. An analytic-numerical method for solving the singular perturbed system of nonlinear PDE related to semiconductor plasma // International Conference "Nonlinear Partial Differential Equations", Alushta, September 15-21, 2003. Book of Abstracts. -P.217.
12. Власов В.И., Безродных С.И. Метод решения сингулярно возмущенной системы нелинейных дифференциальных уравнений // Докл. РАН. - 2004. - 394;№6. - С.731-734.
13. Bezrodnykh S.I. A method for solution of the BVP for singular perturbed system of nonlinear differential equations // International Conference "Differential Equations and Related Topics" dedicated to I.G.Petrovskii. Moscow, May 16-22, 2004. Book of Abstracts. Moscow. - P.28.
14. Безродных С.И., Власов В.И. Метод решения системы нелинейных дифференциальных уравнений, моделирующей полупроводниковый диод // Spectral and Evolution Problems. -2004. - 14. - С.60-72.
15. Безродных С.И., Власов В.И. Краевая задача для моделирования физических полей в полупроводниковом диоде // Ж. вычисл. мат. и матем. физ. - 2004. - 44;№12. - С.2232-2263.
16. Безродных С.И., Власов В.И. Эффективное решение сингулярно возмущенной системы нелинейных дифференциальных уравнений // Современная математика и ее приложения. - 2005. - 36. - С.8-17.
17. Безродных С.И., Власов В.И. Эффективный метод решения сингулярно возмущенной системы нелинейных дифференциальных уравнений // Современная математика. Фундаментальные направления. - 2006. - 15. - С.45-58.
18. Безродных С.И., Власов В.И. Высокоэффективный аналитико-численный метод решения одного класса сингулярно возмущенных систем нелинейных дифференциальных уравнений // Международная конференция по прикладной математике и информатике, посвященная 100-летию со дня рождения академика А.А.Дородницына. Москва, ВЦ РАН, 7-11 декабря 2010 г. Тезисы докладов. - С.87-89.
19. Зи С. Физика полупроводниковых приборов / М.: Мир, 1984.
20. Selberherr S. Analysi and simulation of semiconducter devices / Springer Verlag, 1984.
21. Kano K. Semiconducter devices / Prentice Hall, 1988.
22. Markowich P., Polak S.J., Den Heijer C., Schilders W.H.A. Semiconductor device modeling from the numerical point of view // Internat. J. Numer. Mathods Engng. - 1987. - 24;№4. -P.63-83.
23. Zirkle T.E., Schroder D.K., Backus C.E. The superlinearity of shortcirquit current of low resistance concentrator Solar Cells // IEEE Trans. - 1989. - ED-36;№7. - P.1286-1294.
24. Хэдинг Дж. Введение в метод фазовых интегралов (метод ВКБ) / М.: Мир, 1965.
THE CLASS OF SINGULARLY PERTURBED SYSTEMS OF NONLINEAR DIFFERENTIAL EQUATIONS.
ANALYTIC-NUMERICAL METHOD OF SOLUTION
S.I. Bezrodnykh* **, V.I. Vlasov*
*Dorodnitcyn Computing Centre, Russian Academy of Sciences,
Vavilova St., 40, Moscow, 119333, Russia **Sternberg Astronomical Institute, MSU,
Universitetskii Av., 13, Moscow, 119992, Russia, e-mail: sergeyib@pochta.ru, vlasov@ccas.ru
Abstract. The analytic numerical method is proposed for solving of the certain class of singularly perturbed systems of nonlinear differential equations. The method is based on the operator Newton method, on the parameter continuation method and on the constructing of initial approximation in explicit form. Effectiveness of the method is demonstrated on important special case of systems under consideration which appears when modeling physical fields in semiconductor diode. For this case the Frechet derivative, the Green function of appropriate differential equation are found in analytic form. Numerical realization demonstrated high effectiveness and super-exponential rate of convergence of the proposed method.
Keywords: singularly perturbed systems, analytic-numerical method, physical fields, semiconductor devices.
