CC BY
16
УДК 519.6:533.9:517.9
Д. С. Филиппычев
ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ПОГРАНИЧНОЙ ФУНКЦИИ1
(кафедра автоматизации научных исследований факультета ВМиК)
1. Введение. В работе [1] было рассмотрено применение асимптотического метода пограничных функций (АМПФ) [2-4] к сингулярно возмущенному интегро-дифференциальному уравнению плазма-слой (УПС) [5], описывающему поведение электростатического потенциала и(^) как в основном объеме плазмы, так и в узком пристеночном (ленгмюровском) слое. Здесь безразмерные величины и(^) и £ соответствуют потенциалу и линейной координате. В УПС при старшей (второй) производной стоит малый параметр ц2 (ц <С 1). За вырожденное решение (ц = 0) АМПФ принималось решение интегрального уравнения плазменного приближения щ(^), соответствующее одной из использованных в работе [5] функции источника:
Решение щ(^) получается в результате численного решения уравнения (1) методом деления отрезка пополам. При £ ^ 1/2 оно представляет собой "полочку" постоянной величины щ(^) = const = = и0(1) = щ = 0,40445. Здесь и ниже числовые значения соответствуют выбору параметров задачи работы [5], а также работ [1, 6, 7].
Уравнение, описывающее поведение пограничной функции нулевого порядка (первого члена сингулярной части разложения АМПФ) Y(() = По(С)> гДе С = (1 — — растянутая координата, было выведено в работе [1]:
Здесь с = ехр( — щ), поскольку рассматривается узкая область пристеночного слоя (при V = Ю-2: О ^ £ ^ 10 < 50), в которой ио(С) = их. На решение уравнения (2) накладываются краевые условия У (£ = 0) = Сщ'ч у (С °о) —> 0; Сщ = и\¥ ~ и1] и\¥ = = 1) = 2,96617 — значение решения полного УПС на правой границе [6, 7].
Решение уравнения (2) содержит основную часть падения потенциала в узком пристеночном слое. Поэтому уравнение является фактически уравнением дебаевской экранировки. Вследствие этого (2) и близкие к нему уравнения находят широкое применение при теоретическом исследовании задач физики плазмы (пристеночной, пылевой и т.д.).
В работе [1] для уравнения (2) рассматривалась задача Коши. Результаты расчетов [1] привели к следующему поведению решения: с увеличением значения £ величина У(() сначала, как и ожидалось, уменьшалась, а затем увеличивалась вместо стремления к нулю. В этой работе рассматривается другой путь решения уравнения (2), при котором находится общее аналитическое решение для "модифицированной" первой производной как функции новой независимой переменной. Далее строится разностная схема для численного решения (2). Построенный в работе численный алгоритм решения этого уравнения оказался пригодным как для решения краевой задачи, так и для корректного решения задачи Коши. В последних разделах работы рассматриваются асимптотические решения для пограничной функции и проводится их сравнение с численным решением.
2. Численный алгоритм решения уравнения. Построим алгоритм численного решения уравнения (2). Замена зависимой переменной Z{Cí) = ехр(—У(£)) приводит к уравнению
1 Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект № 02-01-00299.
(1)
(2)
После перехода к новым неизвестным р^) = и = р2 окончательно получаем
г^-2П = -2с(1-г)г\ (з)
Уравнение (3) имеет точное аналитическое решение
11{г) = г2{2с[г-ьп{г)] + с1}. (4)
Из дополнительного условия У (С —> оо) = 0, следует —> оо) = 1, а с учетом экспоненциального
убывания У (С) на бесконечности получаем р = = — . В результате находятся
С—5-00
постоянная С\ = — 2с и решение уравнения (3) с граничным условием У (С —> сю) = 0:
т](г) = 2с2{г- щг) - 1}. (5)
С учетом равенства Ьп(^) = — У (С) решение (5) переписывается в виде = 2с2{^+У — 1}. Значение первой производной решения уравнения (2) вычисляется по формуле
% = -ъ% = -'г1а>=-'гюЛ. (в)
Численное решение уравнения (2) проводилось по следующей схеме, которую в дальнейшем будем называть схемой Э:
1) у =У(0) =
2) ^ = ехр(—У,), Г,3 = 2сг2[гз + У3 - 1];
3) Рз = у/Т], ®з = {<ЖК)з = " ехР(¥з)Рз'>
4) У,+1 = У, + кБ3 + к2[1- ехр(—У,)]с/2, к = 0+1 - С¥Г
Расчеты по приведенной схеме проводились на неравномерной дискретной сетке по £ с минимальным шагом кг$ = 0,015625. Такая сетка ранее использовалась в расчетах [1] и была согласована с частью (10% от правой границы) сетки, рассмотренной в работах [6, 7]. В расчетах использовалось краевое значение Сщ = 2,56172, полученное при решении полного УПС [6, 7]. Результаты расчета представлены в первых колонках табл. 1. В отличие от решения, полученного в [6, 7], результаты расчета по схеме Э показали монотонное убывание функции У (С) Д° значения на правой границе = 10), равном 0,00111. Использование равномерной сетки с шагами к = кг$ , кг$ /2 и кг$ /4 привело практически к тем же самым результатам. Отличия (на единицу) наблюдались только в 5-м значащем знаке.
Таблица 1
УЛз = Сш ■
3 С ПС.,") Аг 6% А* 6% Ащ 6%
2 0,0156 2,5387 2,5292 0,37 2,5422 -0,14 2,5527 -0,55
3 0,0313 2,5158 2,4971 0,74 2,5228 -0,28 2,5436 -1Д
4 0,0469 2,4931 2,4655 1,11 2,5035 -0,42 2,346 -1,7
17 0,2500 2,2110 2,0885 5,54 2,2662 -2,50 2,204 -9,7
33 0,5000 1,8974 1,7027 10,26 2,0048 -5,66 2,869 -20,3
65 1,0000 1,3745 1,1318 17,66 1,5689 -14,14 2,416 -48,4
129 2,0000 0,6812 0,5000 26,60 0,9609 -41,05 1,271 -138
177 3,0000 0,3196 0,2209 30,89 0,5886 -84,11 1,968 -305
209 4,0000 0,1454 0,0976 32,87 0,3604 -147,9 1,335 -610
241 5,0000 0,0651 0,0431 33,77 0,2207 -239,1 0,237 -1165
257 6,0000 0,0289 0,0190 34,19 0,1352 -367,0 0,564 -2168
273 7,0000 0,0128 0,0084 34,39 0,0828 -545,5 0,232 -3979
289 8,0000 0,0057 0,0037 34,49 0,0507 -793,4 0,169 -7246
305 9,0000 0,0025 0,0016 34,55 0,0311 -1137 0,323 -13141
321 10,000 0,0011 0,0007 34,59 0,0190 -1614 0,648 -23771
3. Асимптотические решения рассматриваемого уравнения. Асимптотическое поведение решения уравнения (2) с учетом граничных условий имеет структуру вида
= С№ X ехр(-УАС)- (7)
Результаты вычислений с использованием трех значений параметра А и их сравнение с численным решением уравнения (2) представлены в табл. 1. Приведенные в этой таблице значения параметра А были получены с использованием разложения экспоненты правой части уравнения (2) в ряд Тэйлора по степеням У (С) (параметр А\ = с = ехр( — и\)) и по степеням Сщ — ^(С) (см- [1]) (соответственно Ацг = ехр( — ицг))- Третье значение было получено в работе [10] на основании формализма метода дуального оператора [8, 9] А* = А°(( = 0), где А° (( = Со) = с(1 - ехр(-У(Со))/У(Со)), И0) =
Ниже будут использоваться следующие обозначения:
53 = 100 X (1 — УАв(^)/У(0)) — относительная ошибка в ^'-м узле сетки, где У(С^) — численное решение уравнения (2), полученное по схеме Э;
8л — относительная ошибка на правой границе области (С = 10).
В табл. 1, 2 относительная ошибка обозначена 8%;
Д1 % и Дю% — области, в которых относительная ошибка меньше 1 и 10% соответственно.
Характерные черты рассмотренных асимптотических решений следующие:
Ах = 0,6673: Д1% = (0 < С <~ 0,0625); Д10% = (0 < С < 0,5); 62 = 0,37%; 6Я = 34,59%;
А* = 0,2404: Д1% = (0 < С <~ 0,109); Д10% = (0 < С < 0,78125); 62 = -0,137%; 6Я = -1,6 X 103%;
А\у = 0,0515: Д1% = (0 < С <~ 0,03125); А10% = (0 < С < 0,266); 62 = -0,55%; 6Я = -2,37 X 104%.
Видно, что по величине ошибки во втором узле сетки 62, а также по величине областей % и Длучшей оказывается асимптотика А*. Отрицательное значение относительной ошибки означает, что асимптотическое решение больше численного.
Таблица 2
5ц Локальный максимум д1% дю% ±
а 5% зь Сь 5% 3 С 3 С 3± с±
34,59 4 0,047 33 0,500
а5 25,42 22 4,500 29,36 4 0,047 35 0,531
0,600 0,136 176 2,969 21,54 5 0,063 43 0,656
0,5999 0,071 176 2,969 21,52 5 0,063 42 0,641
0,599 -0,511 175 2,938 21,38 5 0,063 42 0,641 320 9,938
0,590 -6,552 173 2,875 20,05 5 0,063 43 0,656 300 8,688
0,580 -13,75 164 2,594 18,65 5 0,063 46 0,703 284 7,688
0,570 -21,50 157 2,438 17,29 5 0,063 48 0,734 270 6,812
0,500 -96,12 101 1,563 9,13 7 0,094 228 4,594 194 3,531
0,400 -313,7 39 0,594 1,64 16 0,234 150 2,328 80 1,234
0,350 -522,5 11 0,156 0,11 43 0,656 107 1,656 21 0,313
0,300 -865,4 17 0,250 76 1,172
0,250 -1456 9 0,125 54 0,828
а* -1614 8 0,109 51 0,781
0,200 -2538 6 0,078 40 0,609
0,100 -9674 4 0,047 24 0,359
аш -23771 3 0,031 18 0,266
Поведение относительной ошибки 8 асимптотического решения (7) в зависимости от параметра А представлено в табл. 2, в которой показаны:
• относительная ошибка в узле правой границы 8Я = 8(( = (я = 10); номера узлов сетки ^ и их координаты в которых достигается максимальное значение 8^ внутри области расчета, локальный максимум;
• верхняя граница областей Д1% и Дю% (нижняя граница равна нулю) — за верхнюю границу принимается первый узел, в котором ошибка больше 1 и 10% соответственно;
• координата смены знака ошибки ("±") — за ]± принимается первый узел, для которого 82 < 0. Значения параметра А = А^, А* и Ац/ соответствуют значениям, используемым в табл. 1; А5 = А°(С = 5) = 0,646084.
С уменьшением А величина 8Я уменьшается, проходя через нуль при значении равном Ао (А = 0,59979, 8Я = 6,9 X 10"5; А = 0,599789, 8Я = -5,77 X 10"4). После этого отрицательные значения возрастают по абсолютной величине до больших значений. Например, при А = Ац/ = 0,0515 — 5Я = -23771%.
Использование 8Я в качестве величины, характеризующей близость асимптотического и точного решений, оправдано только в том случае, когда 8Я является монотонной функцией координаты, т.е.
внутри области Sj монотонно возрастает по абсолютной величине. Однако, по данным табл. 2, имеется область параметра A G [0,6 -т- 0,58], в которой ошибка максимальна внутри области Sl > Sr. Область немонотонности начинает появляться при А = Alb = 0,6665, когда = 320 является предпоследним узлом сетки, a Sl = Sr = 34,25%. По мере уменьшения А локальный максимум ошибки смещается к началу координат, а его величина уменьшается, оставаясь всегда положительной. При этом Sl > Sr до достижения параметром А значения Arl, при котором величина Sr, возрастая по абсолютной величине, становится равной Sl = После этого доминирующей ошибкой стано-
вится SR. Значение Акь и 0,574554. При этом SR = -17,9011%; SL = 17,9006%. Для А = 0,574555 соответственно Sr = —17,9004%; Sl = 17,9007%. В обоих случаях = 161; (l = 2,5. Заметим, что при Arl достигается лучшее приближение асимптотики к точному решению: ошибка по всей области расчета составляет приблизительно 18%.
Отрицательные значения Sr появляются при А ^ Ао- Изменение знака S происходит в последнем (граничном) узле сетки j± = J321. По мере уменьшения А узел j± смещается влево. После того как j± достигает второго узла сетки при А = Ale (^4 = 0,3349991; j± = 3, S2 = 0,00113%; А = 0,334999; j± = 2, ¿2 = 0,000893%), внутри области расчета ошибка монотонно увеличивается по абсолютной величине, достигая максимального значения на правой границе области расчета Sr.
Верхние границы областей, в которых ошибка меньше 1 и 10%, представлены в колонках табл. 2, обозначенных соответственно Ai% и Дю%- В расчетах были получены следующие длины областей: maxAi% = 0,65625 для А = 0,35; тахА10% = 4,594 для А = 0,5. Как отмечалось выше, когда параметр А изменяется в интервале {Ale, ^4о), внутри области расчета относительная ошибка S проходит нулевое значение. Поэтому в окрестности точки S = 0 могут возникнуть дополнительные области Ai % и Аю%- Длина дополнительной области Ai% не превышает величины 0,5. Следует отметить, что при приближении j± к началу координат дополнительные области могут и не появиться. Для этого необходимо и достаточно, чтобы локальный максимум S не превышал уровня соответственно в 1 и 10%, т.е. когда соответствующая область содержит точку изменения знака j±. Иллюстрацией к сказанному может служить строка табл. 2 для А = 0,35.
4. Численное решение уравнения при выборе постоянной С\ для задачи Коши. В
общем случае в схеме S вычисление r¡j осуществляется по формуле (4), которую с использованием соотношения Ln (Z) = — У (С) можно переписать как r¡j = {2 c(Zj + Yj) + C\)Z2-. Значение постоянной С\ вычисляется в результате рассмотрения приведенного выражения в краевой точке (о = 0:
С\ = щ/Z2 - 2c(Z0 + Yo) = {dY/dQl - 2c(Z0 + Y0). (8)
При постановке задачи Коши для уравнения (2) задается величина первой производной в начале координат Do = (dY/dQ£=q. После подстановки в (8) выражения щ — = ^о^о föo = ехР(—Cw)), вытекающего из (6), определяется значение постоянной С\ \ С\ = Dq — 2c{Zq + Cw)-
Расчеты по схеме S проводились для двух значений Dq\
1) D0 = -fiQ = -0,9833 [1];
2) Dq = —Cw"*/~Ä= —1,2560, А = c( 1 — exp(—Cw))/Cw — выражение, полученное в работе [10] на основании формализма дуального оператора.
Результаты расчетов сравнивались с "точным" решением, за которое была выбрана разность между численным решением УПС [6, 7] и решением уравнения плазменного приближения вблизи правой границы и(() = щ. Для обоих рассмотренных случаев полученные значения вблизи начала координат меньше "точного" решения (ошибка положительна). В дальнейшем в точке (± происходит смена знака ошибки, что означает более медленное убывание численного решения по сравнению с "точным". После смены знака ошибка монотонно возрастает по абсолютной величине. Начиная с точки (с решение практически выходит на постоянное значение Ye, с чем и связаны большие значения относительной ошибки. Максимальная ошибка ¿max достигается на правой границе области расчета. Были получены следующие результаты:
• для варианта 1: S2 = 2 х 10"3%; С± = 0,0469; Сс = 1,672; Yc = 1,739; ¿max = 2,7 X 103%; Ai% = [0,0,306]; Аю% = [0,1,031];
• для варианта 2: S2 = 0,17%; С± = 1,8125; (с = 2,359; Yc = 1,136; ¿max = -1,7 X 103%; Ai% = [0, 0,109; 1,750,1,859]; Д10% = [0, 2,156].
Видно, что при различных значениях производной решение стремится к различным постоянным. Оценивая полученные результаты, надо иметь в виду, что была рассмотрена только пограничная
функция нулевого порядка. Большим значениям относительной ошибки соответствуют относительно небольшие значения абсолютной ошибки (порядка 0,1). Такое различие в значительной мере может быть обусловлено пограничной функцией первого порядка, поскольку в работе [6, 7] использовалось значение ¡i = Ю-2.
5. Заключение. В данной работе построен эффективный алгоритм численного решения дифференциального уравнения второго порядка, описывающего поведение пограничного слоя (2). Численный алгоритм является корректным как при решении задачи Коши, так и при задании второго условия на бесконечности. В первом случае в зависимости от значения заданной производной решение выходит на различные постоянные значения. Кроме того, в работе представлены асимптотические решения уравнения (2), которые достаточно хорошо описывают поведение решения вблизи начала координат.
СПИСОК ЛИТЕРАТУРЫ
1. Филиппычев Д. С. Метод пограничных функций для получения асимптотического решения уравнения плазма-слой // Прикладная математика и информатика. № 19 / Под ред. Д.П. Костомарова, В.И. Дмитриева. М.: МАКС Пресс, 2004. С. 21-40.
2. Васильева A.B., Бутузов В.Ф. Асимптотические разложения решений сингулярно возмущенных уравнений. М.: Наука, 1973.
3. Васильева А. Б., Бутузов В.Ф. Сингулярно возмущенные уравнения в критических случаях. М.: Изд-во МГУ, 1978.
4. Васильева А. Б., Бутузов В. Ф. Асимптотические методы в теории сингулярных возмущений. М.: Высшая школа, 1990.
5. Emmert G.A., Wieland R.M., Mense А.Т., Davidson J.N. Electric sheath and presheath in a collisionless, finite ion temperature plasma // Phys. Fluids. 1980. 23. N 4. P. 803-812.
6. Филиппычев Д.С. Численное моделирование уравнения плазма-слой с использованием сгущающейся сетки // Прикладная математика и информатика № 14 / Под ред. Д.П. Костомарова, В.И. Дмитриева. М.: МАКС Пресс, 2003. С. 35-54.
7. Филиппычев Д.С. Численное моделирование уравнения плазма-слой // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2004. № 4. С. 32-39.
8. Cacuci D.G., Perez R. В., Protopo pese u V. Duals and propagators: A canonical formalism for nonlinear equations // J. Math. Phys. 1988. 29. N 2. P. 335-361.
9. Cacuci D.G., Karakashian O.A. Benchmarking the propagator method for nonlinear systems: A Burgers-Korteweg-de Vries equation //J. Comput. Phys. 1990. 89. N 1. P. 63-79.
10. Филиппычев Д.С. Применение формализма дуального оператора для получения пограничной функции нулевого порядка уравнения плазма-слой // Прикладная математика и информатика. № 22 / Под ред. Д.П. Костомарова, В.И. Дмитриева. М.: МАКС Пресс, 2005. С. 76-90.
Поступила в редакцию 07.06.05
