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

УДК 517.583

О ВЫЧИСЛЕНИИ ЭЛЛИПТИЧЕСКИХ ИНТЕГРАЛОВ ТРЕТЬЕГО РОДА

С.Д. Симонженков

The algorithm for numerical evaluating the Legendre’s incomplete integral of the third kind with parameter (character) less than one is described. The Fourier cosinus-expansion of the integrand part involving a character is used. Termwise integration gives a series for which summation a Clenshaw’s technique is applicable

1. Введение

При решении некоторых задач физики, техники, геодезии и др. возникает необходимость вычисления эллиптических интегралов, в частности, интегралов третьего рода

П = П(<д, h,k) = J(1 — k2sin2t)~ll2 (1 — hsin2t)~1dt, о

параметры которых удовлетворяют условиям

0 < < тг/2, 0 < к < 1, h < 1, h ф- 0. (1)

Среди соответствующих многочисленных вычислительных методов укажем следующий, основанный на разложении подынтегрального множителя fit) = (1 — к2sin2t)~ll2 в ряд

ОО

/«) = £'<■ ncos2nt. (2)

п=О

(Здесь и далее штрих в знаке суммы означает, что при п = 0 слагаемое берется с коэффициентом 0.5). Коэффициенты ап в (2) задаются через гипергеометрическую функцию Гаусса в виде [1, с. 34]:

ап = 2 (-!)"(! + q) qn (l/2)n 2F, (n + 1/2; 1/2; n + 1; q2)/n !

0 2000 С.Д. Симонженков

Омский государственный педагогический университет

и удовлетворяют соотношениям

(п + l/2)an+i + 2n(2/k2 - 1 )an + (n - l/2)an_1 =0, n > 1, = 1. (4)

n=0

В формуле (3) q = (1 — л/1 — k2)/( 1 + л/l — &2), а выражение (...)„ обозначает символ Похгаммера. Почленным интегрированием из (2) можно получить разложение в ряд по синусам четных дуг эллиптического интеграла первого рода

F = /(1 — k2sin2t)~ll2dt. Аналогичным образом получается соответствующее о

разложение для интеграла второго рода Е = J(1 — k2sin2t)l^2dt.

о

Применительно к интегралу П разложение (2) используется следующим образом. Из (2) получаем, что

ОО

п = У^/anPn,

п=0

(5)

ч>

где pn = f (1 — hsin2t)~1cos2nt dt. Поэтому вычисление П сводится к нахожде-о

нию величины Пдг = ^}2n=o ' anPn при достаточно большом N. Выбор N зависит от задаваемой погрешности г и, как нетрудно проверить, может быть осуществлен исходя из условия 2p0qN+1 / (1 — q)2 < е. С другой стороны, интегралы рп удовлетворяют рекуррентному соотношению

Рп+1 + anpn + finpn-1 = n > 1, (6)

где an = 4/h — 2, (3n = 1, yn = (2/h)(sin2rup)/n, n > 1. Поэтому вычисление Пдг можно осуществить суммированием Кленттто [1, гл. 11] по следующему алгоритму. Задаем Ндг+i = 0, Ндг+2 = 0 и счетом назад находим

Вп оспВп+\ [3n+iВп+2 -)- £пап^ п IV, N 1,..., 0, (7)

где еп = 1 при п > 0 и е0 = 1/2. После этого окажется, что

Пдг = (Во + a0Bi)p0 + PiBi + сгдг,

где сгдг = Yln=o Bn+iln, если доопределить у0 = 0. Здесь последовательность {in} удовлетворяет соотношению

7п+1 + а/7n + P'nln-i = о, n > 1, (8)

с коэффициентами a'n = — (2n cos2p)/(n + 1), fi'n = (n — 1)/(га + 1), n > 1. Следовательно, сгдг вычисляется аналогично: при Адг+i = Адг+2 = 0 надо найти

Ап =—anAn+i — l3'n+1An+2Bn+i, п = N, N — 1,..., 0. (9)

Отсюда получаем, что сгдг = (+о + сф+Ойо + А71 = А71. Окончательно,

Пдг — (В0 + a0Bi)po + Bipi + А171, (10)

где

Ро

л/1 — h

arctg(\/l — htgp), щ

2 2 2 С1 ~~ л)№ + hip' 71 = hsm 2lp'

Таким образом, вычисление П сводится к нахождению элементарной функции ро - в этом достоинство описанного подхода. Его недостатком является тот факт, что коэффициенты а0, од,..., адг подлежат нахождению. Например, на основе (4), их можно вычислить устойчивой обратной рекурсией с помощью алгоритма Миллера [1, гл. 12].

2. Формулировка результата

В настоящей работе описанный метод используется с той разницей, что в ряд (2) разлагается другой подынтегральный множитель f(t) = (1 — hsin2t)~l. На этот раз коэффициенты разложения находятся явно, т.е. an = 2( — l)nqn/\/l — й, где q = (1 — \/1 — h)/(1 + \/1 — й), а в формуле (5)

pn = J(1 — k2sin2t)~1!2cos2nt dt. (11)

о

Как и выше, ряд (2) сходится не медленнее геометрической прогрессии со знаменателем д, при этом

П

Пдг| <

2ро

Vl - h 1 - \q\

ы

JV+1

В равенстве (6)

2 п

an = 2 (— — 1)--------—

Ук2 п + 1/2

, а = AV % = А а - ^»а>1/2- («I

п + 1/2’ к2 п + 1/2

Коэффициенты удовлетворяют (8) при следующих значениях аф, /3/:

. га+ 1/2 „ п. п- 1/2

V = —2--------соз2(д, (Зп =

п

+ 3/2

п

+ 3/2

Оказывается, что в (10)

Ро = F, pi = (1 — 2/k2)F + (2/к2)Е, 71 = (4/ЗА;2)7/1 — к2sin2sin2<p. (13)

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

вычисления существуют компактные и эффективные алгоритмы [2, гл. 17]. Поэтому предполагается, что в излагаемой ниже схеме интегралы F, Е считаются известными.

Перейдем к описанию вычислительной схемы. Для нахождения П с погрешностью е при аргументах, удовлетворяющих условиям (1), надлежит сделать следующее (используются обозначения предыдущего пункта).

1. Вычислить q, подобрать целое N > 0 такое, чтобы

2 F

Vl - h 1 - \q\

< £.

2. Сформировать массивы {_£?„}, {Ап}, п = 0,1,..., А + 2 в соответствии с формулами (7), (9).

3. Вычислить Пдг, используя (10), (13).

Для иллюстрации приведем следующий пример. При е = 10-8 найти

П(30°, 0.5, sin 60°). Согласно [2] (табл. 17.5, 17.6), имеем, что F = 0.54222911, Е = 0.50609207. Имеем далее, что q = 0.17157288, N = 10, aw = 2.2... • 10“8. Результат П = 0.56836557 согласуется с табличным значением [2] (табл. 17.9).

Замечание . Применять предложенный алгоритм нельзя, если параметр к близок к нулю либо h близко к 1 или h < 0 и велико по модулю. Действительно, при к ss 0 величина 2 (2/к2 — 1) велика; в этом случае Пдг находится как разность больших величин, что приводит к потере значащих цифр. Например, при вычислении П(90 °, 0.5, sin 15 °) с погрешностью е = 10-5 надо брать N = 6. При вычислениях на 8-разрядном калькуляторе имеем 2(2/к2 — 1) = 57.71282, В\ = -19181.611, В0 = -164.81929, Pl = -1.38504 • 10“2, р0 = 1.598142. Получим П = роВ0 -\- pi Вi = 2.26836, что плохо согласуется с табличным значением П = 2.26685. При h -У —оо или h —У 1 имеем, что |д| —у 1. Следовательно, в алгоритме N велико. Его применение наиболее эффективно, если |д| < 1/2, что соответствует области —8 < h < 8/9. Как показывают вычисления, при таких h и к2 > 1/2 алгоритм «работает» удовлетворительно. Существенно используемое разложение (2) для fit) = (1 — hsin2t)~l можно найти в руководствах по рядам Чебышева (см., например, [3, с. 144]). Вывод равенства (6) для интегралов (11) с коэффициентами (12) имеется в [4].

Литература

1. Люк Ю. Специальные математические функции и их аппроксимации. М.: Мир, 1980.

2. Справочник по специальным функциям / Под ред. М. Абрамовича и И. Стиган.

М.: Наука, 1979.

3. Пашковский С. Вычислительные применения многочленов и рядов Чебышева. М.: Наука, 1983.

4. Литвин А.И., Симонженков С.Д. О вычислении некоторых эллиптических интегралов. Деи. в ВИНИТИ Омским с.х. институтом. 12.12.90. N 6202-В 90.