ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 5, с. 798-806
УДК 519.624.2
ВЫЧИСЛЕНИЕ СФЕРОИДАЛЬНЫХ ФУНКЦИЙ I РОДА ДЛЯ КОМПЛЕКСНЫХ ЗНАЧЕНИЙ АРГУМЕНТА И ПАРАМЕТРОВ
© 2015 г. А. А. Абрамов*, **, Е. Д. Калинин*, **, С. В. Курочкин*
(* 119991 Москва, ул. Вавилова, 40, ВЦ РАН; ** 141700Долгопрудный, М. о., Институтский пер., 9, МФТИ) e-mail: alalabr@ccas.ru; e.kalinin@inbox.ru; kuroch@ccas.ru Поступила в редакцию 06.11.2014 г.
Предложены методы: а) нахождения собственных значений волнового сфероидального уравнения с комплекснозначными параметрами, расположенных в заданной области комплексной плоскости; б) вычисления значений соответствующей функции для комплексных значений аргумента. Библ. 8. Фиг. 2. Табл. 2.
Ключевые слова: волновые сфероидальные функции, спектральная задача, рекуррентные формулы, численная устойчивость.
DOI: 10.7868/S0044466915050038
ВВЕДЕНИЕ
Сфероидальные функции играют большую роль в математической физике (см., например, [1]). В настоящее время существенной является задача вычисления этих функций для комплексных значений параметров и аргумента (см., например, [2]—[5]). Для этого случая в настоящей работе рассматриваются следующие вопросы: вычисление собственных значений возникающей спектральной задачи и вычисление сфероидальной функции I рода при значениях аргумента, не превосходящих или несущественно превосходящих по абсолютной величине единицу.
1. ПРЕДВАРИТЕЛЬНОЕ ИССЛЕДОВАНИЕ
В комплексной области рассматривается уравнение
j(( 1 - z ) j) + (х + 40( 1 - z ) - и = 0, (1.1)
dZ dZ 1_z
где X, 0, m — комплексные параметры уравнения, 0 и m заданы. Это уравнение имеет три особые точки: +1 и —1 — регулярные точки с показателями ±m/2, и да — иррегулярную точку. В регулярных точках соответствующие два линейно независимых решения ведут себя следующим образом:
и1,2(z) ~ (1 _ z2)±m/2 при m Ф 0 ,
u1(z) ~ const, и2(z)~ ln(1 _ z ) при m = 0.
Так как в уравнение входит только m2, то без ограничения общности можно считать, что Re(m) > 0. Далее предполагаем, что Re(m) > 0 или m = 0. Случай Re(m) = 0, m Ф 0, не рассматривается, так как тогда каждое решение уравнения (1.1) ограничено в окрестности точек +1 и —1 и рассматриваемая далее спектральная задача не возникает.
Сфероидальная функция I рода определяется как нетривиальное решение уравнения (1.1), ограниченное в окрестности точек + 1 и —1. Эта функция может быть неоднозначной. При заданных 0 и m тем самым возникает спектральная задача, в которой X — спектральный параметр.
Переходя к новой искомой функции (см. [1])
v( z) = и (z) / (1 _ z2 )m/2, (1.2)
получаем уравнение
(1 - г2К' - 2(т + 1)zV + [X - т2 - т + 40( 1 - г2)] V = 0. (1.3)
Требование ограниченности функции ы^) в окрестности точек +1 и —1 эквивалентно тому, что v(z) — целая функция.
Далее будут использоваться функции Рт + г (¿), г — целое, связанные с присоединенными функциями Лежандра (см. [1], [6]) Рт + г следующим образом:
Рт + г (г ) —
2г"г2(т+11))рт+г(г), г — о, 1,2,...,
Г( 2 т + 1) 0, г — -1, -2,..
В частности,
Рт(г) — (1 - г2)т/2, Рт + 1(г) — (2т + 1)(1 - г2Г\
Из известного рекуррентного соотношения для Р^ + г и формул для Рт и Р„, + г (z) следует, что функции Рт + г (¿) удовлетворяютрекуррентномусоотношению
(г + 1)Рт+г +1 (г) — (2т + 2г + 1)гРт + г(г) - (2т + г)Рт+г - 1(г) для целых г, г = 0, ±1, ±2, ... .
Далее, в соответствии с [1], функция ы^) будет представляться разложением в ряд по функциям рт + г (¿), а v(z) — по соответствующим многочленам рт + г где Рт + г (z) = (1 - z)2pm + г (z) (в частности, рт (z) = 1, рт + г = (2т + 1)z), для которых рекуррентная формула, очевидно, та же:
(г + 1 )рт + г +1 (г) — (2т + 2г + 1)грт + г(г) - (2т + г)рт + г-1 (г). (1.4)
Такое разложение имеет место и в том случае, когда т, 0, X, z комплексные.
Уравнения (1.1) и (1.3) не меняются при замене z на —z, поэтому достаточно рассматривать случаи:
и (г), v(z) — четные функции, (1.5)
и (г), v(г) — нечетные функции. (1.6)
Тогда получаем (см. [1])
V (г) — ^ +г( г), (1.7)
где суммирование ведется по г = 0, 2, 4, ... для случая (1.5) и по г = 1, 3, 5, ... для случая (1.6). Коэффициенты разложения удовлетворяют следующей рекуррентной формуле (для обоих случаев):
ла+2 + (вг - хк + са - 2 — о, (1.8)
где
л — 4 0 ( 2 т + г + 2 ) ( 2 т + г + 1 ) г (2 т + 2г + 3 ) ( 2 т + 2 г + 5 ),
Вг — (т + г)(т + г + 1) - 8 0((г + т) (г + т + 1 ) + т2 - 1 ), г > 1, г У ' (2т + 2г + 3)(2т + 2г - 1)
о ( , 1) 80(т+1)
Во — т(т + 1) -—5.--1
2т + 3
(1.9)
(упрощение (2т2 + т — 1)/(2т — 1) = т + 1 существенно при т = 1/2),
= -4Мг-!)-, г > 2, с0 = С1 = 0.
г (2т + 2г - 3)(2т + 2г - 1) 0 1
При г = 0 и г = 1 формула (1.8) принимает вид
ЛЛ+2 + В - X) йг = о. (1.10)
При 0 Ф 0, если не учитывать (1.10), то в двумерном пространстве возникающих для (1.5) и (1.6) последовательностей йг некоторое одномерное пространство содержит последовательности, стремящиеся к нулю очень быстро: г24г+ 2/4г —»- — 0; остальные последовательности очень быстро растут: (4г+1/4г)/т2 —»- — 0. Собственные значения X выделяются совокупностью условий: (1.10) и ИШг ^ ^г+2/4 = 0.
Если нужное собственное значение вычислено, то для вычисления v(z) (и затем по формулам (1.2) и(г)) используется ряд (1.7). Его члены вычисляются с использованием рекуррентных формул (1.4), (1.8) (важные детали такого использования см. ниже в разд. 3), частичные суммы этого ряда — многочлены.
При больших г рекуррентная формула (1.4) имеет предельный вид
Рт + г + 1( = 2%Рт + г ) - Рт + г - 1 ),
откуда следует, что для больших \г\ значения рт + г (г) растут не быстрее \2г\г. Поэтому ряд (1.7) равномерно сходится в любой ограниченной области комплексной плоскости и тем самым действительно дает целую функцию.
Отметим, что если 0 = 0, то мы получаем ответ без всех описанных вычислений: для любого целого неотрицательного г получаем X = (т + г)(т + г + 1), у(г) = рт + г (г).
Далее полагаем, что 0 Ф 0.
Для представления искомой функции v(z) мы использовали ряд (1.7). В настоящее время используются более общие разложения (см., в частности, [3], [4], где разложение содержит свободный параметр V), ряд (1.7) и остальные используемые нами формулы — частный случай, получаемый, если взять V = т и положить рт + г (г) = 0 при целом отрицательном г. При этом некоторые формулы упрощаются, в частности, не возникает рядов, в которых индекс суммирования меняется от —да до +да.
2. ВЫЧИСЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИИ
Представим 4г = - 0г/2 с1г для случая (1.5), 4Г = 1 0(г -1)/2 Сг для случая (1.6). Мы получим г! г!
ЯЯ+2 + В4г + СгЯ - 2 = о, (2.1)
где при г > 2 имеем
Л = _Л г0_ В = В г - Х С = С-
г = ( г - 1) г(г + 1) (г + 2 ), г = г( г - 1), г =0 ,
Яо = , Во = В0 - X, Л1 = —, В1 = В1 - X, Со = С = 0. 2 6
Здесь Аг —► 0, Br —► 1, Cr —1 при r—► да. Нужно выделить последовательность, для которой
lim = -1.
Для каждого из случаев (1.5), (1.6) фиксируем гш, одно и то же для всех исследуемых далее X, четное для случая (1.5), нечетное для случая (1.6). Положим
С1гт = 1, Сг„ - 2 = -1.
г^ +» d
г
По формуле (2.1) вычисляем -4, -6, .... Эти вычисления проводятся численно устойчиво, так как последовательности, не пропорциональные той, которую мы хотим вычислить, для больших г растут слева направо и, следовательно, справа налево убывают; численная погрешность в
выборе drm и в выкладках убывает для больших г при счете справа налево. Вычисления реализуемы, так как Сг Ф 0 для всех г > 2.
Для взятого Х окончательно получим
W(X) = + Шо для (1.5), ЩХ) = Лк^ + для (1.6).
Очевидно, Ж(Х) в обоих случаях (2.2) — аналитическая функция. Соотношение Ж(Х) = 0 предельно при гш —»- да эквивалентно тому, что X — собственное значение.
В случае если 0 вещественное и т целое, собственные значения задачи являются вещественными; они естественно нумеруются (см. [1]) и естественной является задача: найти собственное значение с заданным номером (метод решения этой задачи см., например, в [7]). В работе [5] предложен метод вычисления собственных значений Х(0) в случае целых неотрицательных т и комплексного 0. В рассматриваемом сейчас общем случае нет способа как-либо "занумеровать" собственные значения. Поэтому задача вычисления собственного значения с заданным "номером" заменяется следующей задачей. Зафиксируем в комплексной плоскости ограниченную область Q с кусочно-гладкой границей Г. Предположим, что Г не содержит собственных значений (что будет иметь место в "общем положении"). Нужно найти все собственные значения, лежащие в Q.
Предлагаемый метод решения этой задачи состоит в следующем. Имея возможность вычислить значения функции Ж(Х) на Г и пользуясь известными формулами теории функций комплексного переменного, составим систему уравнений для собственных значений, лежащих в Q.
Так как Ж(Х) не обращается в 0 на Г, то имеют место следующие формулы (см., например, [6]):
1 р ащ W(X) = п, (2.3)
г
п
2П |х^1и W(X) = V Х^. (2.4)
г
г = 1
Здесь п — число собственных значений (т.е. нулей функции ЩХ)), лежащих в Q; Х1, ..., Хп — эти собственные значения (кратные нули функции Ж(Х) считаются нужное число раз), к = 1, 2, ...; мы берем далее к = 1, 2, ... п после вычисления п; при интегрировании контур Г проходится так, что Q остается слева.
Далее рассматриваем случай, когда Г состоит из одной замкнутой кривой; обобщение на случай, когда Г состоит из нескольких таких кривых, очевидно.
На Г выбираем сетку р0, ... = р0. Заменив интеграл (2.3) соответствующей суммой, вычислим число ¿0,:
^ = ± V ащ ^) . (2.5)
0 1) ( )
s = 1
Точки р0, ... = для однозначного вычисления -1) должны для всех 5 удо-
влетворять условию |аг§Ж(р')/< п при р' и р", лежащих на дуге . Если нет собствен-
ных значений, близких к Г, то сетку можно взять достаточно грубой. Мы получим п = ¿0. После этого для к = 1, ..., п, аппроксимируя (2.4), вычисляем
М к к ли ч
^ = ^ V Ьк±Нк_11п. (2.6)
к 2 W(ц,_ 1)
5 = 1
Здесь мелкость сетки определяется уже не только требованиями, которые указаны при вычислении Б0, но и требованием получить необходимую малую погрешность в ответе. Отметим высокий порядок то
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.