О табулировании с высокой точностью нулей функций Бесселя

Известия Тульского государственного университета Естественные науки. 2013. Вып. 1. С. 132-141

= ИНФОРМАТИКА =

УДК 519.632.4

О табулировании с высокой точностью нулей функций Бесселя

С. Д. Алгазин

Аннотация. Приводится алгоритм вычисления нулей функций Бесселя, который позволяет вычислить несколько первых десятков нулей с 15-30 знаками после запятой.

Ключевые слова: уравнение Бесселя, численные алгоритмы без насыщения.

Введение. С появлением современных ЭВМ задача табулирования изменилась. Примером нового подхода к задаче табулирования служит справочник Ю. Люка (Yudel L. Luke) [1]. В предисловии К.И. Бабенко к этой книге сказано: «... Раньше, во времена ручных вычислений, мы были, как правило, заинтересованы в том, чтобы иметь обширную таблицу, которая позволяла бы вычислять значение функции с помощью простейшей интерполяции. Теперь, используя при вычислениях ЭВМ, мы очень часто заинтересованы в том, чтобы иметь таблицу минимального объёма, даже если алгоритм восстановления функции достаточно сложен. Таким образом, возникает задача построения оптимальных по объёму таблиц». Далее показано, что возможно вычислить несколько десятков первых нулей функций Бесселя с 30 знаками после запятой.

1. О табулировании нулей функции Бесселя J0

Рассматривается задача на собственные значения для уравнения Бесселя нулевого индекса:

У

Краевая задача (1)—(3) эквивалентна интегральному уравнению

4—) = 2j Ч—-^г) — Ч ~)

-1

G(x,0 = - ln[(x + с + |£ - x|)/2].

Применим для функции [(£+1)/2]у[(£+1)/2] интерполяционную формулу вида

T

(Pnf )(x) = V f (Tk )lk (т) + Rn(x; f), lk (т) = 7------—Г, k = 1,2,...,n,

k=

Tn(x) = cos(n arccos т), Tk = cos[(2k — 1)^/2n], где Rn(x; f) — погрешность интерполяции. Проводя вычисления, получаем

n

Vj = ^Bjk Vk + rn(Tj; y), Bjk = Bk (Tj), yk = y(Tk), k,j = 1,2,...,n, k=i

где

Bk(т) = k+i j, £±i) lk(c)*,

-1

n.(T,y) = I ja ( ^, Ц2-1) Rn (; (i±±)V) >k.

-1

Отбрасывая в (5) погрешность дискретизации, получаем приближенную задачу на собственные значения y = \By. Здесь у — вектор приближенных значений искомой собственной функции у(х) в узлах сетки, Л — приближенное собственное значение. Легко видеть (см., например [2, с. 189]), что max |rn(x; y)| ^ c |Л| (1 + wn)En(y), где с — абсолютная постоянная,

(ши = 0(ln(n)) — константа Лебега интерполяции, а En(y) — наилучшее приближение функции у многочленом степени не выше (п — 1) в норме С). Заметим далее, что собственные функции задачи (1)—(3) целые, а поэтому lim ПEn(y) = 0 (см. [3, с. 254]), т.е. величина погрешности дискретизации

n—

очень быстро стремится к нулю. Возмущение, вносимое в собственные значения отбрасыванием погрешности дискретизации, будет оценено ниже. Здесь обсудим результаты численных расчетов. При п = 5 получим первое собственное значение с четырьмя знаками после запятой, а третье собственное значение - с одним знаком после запятой. При n = 20 первое собственное значение вычисляется с 22 знаками после запятой, а 14-е собственное значение вычисляется с одним знаком после запятой. Последний расчет проводился на ЭВМ БЭСМ-6 с использованием двойной точности (длина мантиссы 80 бит). Время расчета 4 мин. 40 с. вместе с вычислением матрицы. В [4] опубликована программа BESSEL, по которой проводились эти расчеты, а также результаты расчета собственных функций краевой задачи (1)—(3).

Эти результаты получены автором в соавторстве с К.И. Бабенко давно, целью настоящей работы является их уточнение с использованием совре-

менных вычислительных средств. Сейчас появился транслятор с Фортрана (Intel Visual Fortran 9.1), который позволяет вести расчёты с учетверенной точностью REAL*16. С использованием этого транслятора были проведены расчёты с числом точек N = 3-23, N = 110. Эти расчёты сравнивались с таблицей VI работы [5]. Ниже приведены значения \/Л^, i = 1, 2,..., 13 на сетке из N = 23 узлов:

N = 23; EPS= 0.33E-28

2.40482555769577276862163187932315 5.52007811028631064959995531048766 8.65372791291101229166741865684604 11.7915344390140843686267083669051 14.9309177084599940892731158995010 18.0710639711203998901712950375023 21.2116367066011374633859922408520 24.3524696616047016958823129522958 27.4934601396705964187069258148688 30.6347976897899335897734472259312 33.7759497154294753471218322997414 36.9231498885250203330212676349967 39.9510270188728602686403023637232

Курсивом выделены знаки, совпавшие с таблицей VI работы [5]. Таким образом, сетки из N = 23 узлов достаточно, чтобы получить первое собственное значение со всеми знаками из указанной выше таблицы (29 знаков после запятой). Расчёты на сетке N = 110 дают 32 собственных значения (всю первую колонку таблицы VI) со всеми знаками таблицы (28-29 знаков после запятой).

Далее проводились расчёты на сетках из N = 3-23 узлов. Ниже приведена погрешность определения первого собственного значения: 3) 0.11;

4) 0.53e-02; 5) 0.37e-03; 6) 0.54e-04; 7) 0.2e-05; 8) 0.23e-06; 9) 0.64e-08; 10)

0.68e-09; 11) 0.15e-10; 12) 0.14 e-11; 13) 0.27e-13; 14) 0.23e-14; 15) 0.36e-16; 16) 0.28e-17; 17)0.39e-19; 18) 0.28e-20; 19) 0.35e-22; 20) 0.23e-23; 21) 0.26e-25 ; 22) 0.15e-26; 23) 0.33e-28. По этой таблице подбиралась аналитическая зависимость е = e(N). Получена зависимость lne=a+bN+ cN2+dN3+eN4+fN5; a=0.39307047, b=0.58615539, c=0.40323914, d=-0.599229377, e=0.12421718, f=-0.0076898089. Таким образом, скорость убывания погрешности экспоненциальная. Заметим, что метод конечных разностей и метод конечных элементов дают только степенное убывание погрешности.

2. Вычисление функций Бесселя целого индекса

Уравнение Бесселя имеет вид:

1 /к\2

—(v"(r) +—v'(r)) + ( -j v(r) = Лу(т), v(1) =0, |v(0)| < то.

Поэтому легко видеть, что матрица дискретной задачи (матрица, у которой необходимо вычислить собственные значения) имеет вид:

Ло = В-1,

где В — матрица дискретного оператора (обратного к оператору Бесселя), построенная в [4];

Лд = Л0 + к2 К, к = 1,2,... ,т;

К — диагональная матрица с числами (1/г^)2ти = (1 + еоз((2и — 1)п/2т))/2,

V = 1, 2,... ,т на диагонали.

Рассмотрим результаты численных расчётов. На сетке из 23 узлов получено:

МАСНЕР = 1.92Е-4; т = 23, к=1; уХ собственные значения:

1) 3.83170597020751231561443588605205

2) 7.01558666981561875353120866886788

3) 10.1734681350627220491724007823358

4) 13.3236919363142570433728670520907

5) 16.4706300508824154776364636843869

6) 19.6158585101973123876830625220157

7) 22.7600843701865157377758772901100

8) 25.9036722639552867175841920011564

9) 29.0468258315682780761054280922193

10) 32.1898292345730302050572370405816

На сетке из 140 и 280 узлов получено: т = 140 к = 1; у/\Ї собственные значения

1) 3.83170597020751231561443588630900

2) 7.01558666981561875353704998147703

3) 10.1734681350627220771857117767758

4) 13.3236919363142230323936841269486

5) 16.4706300508776328125524604709891

6) 19.6158585104682420211250658841377

7) 22.7600843805927718980530051521825

8) 25.9036720876183826254958554459800

9) 29.0468285349168550666478198835323

10) 32.1896799109744036266229841044603

11) 35.3323075500838651026344790225194

12) 38.4747662347716151120521975577165

13) 41.6170942128144508858635168050605

14) 44.7593189976528217327793527132122

15) 47.9014608871854471212740087225073

16) 51.0435351835715094687330346332243

17) 54.1855536410613205320999662145337

18) 57.3275254379010107450905042437506

19) 60.4694578453474915593987498083826

20) 63.6113566984812326310397624178737

21) 66.7532267340984934153052597500424

22) 69.8950718374957739697305364354993

23) 73.0368952255738348265061175690918

24) 76.1786995846414575728526146235348

25) 79.3204871754762993911844848724880

26) 82.4622599143735564539866106487815

27) 85.6040194363502309659494254933808

28) 88.7457671449263069037359164348545

29) 91.8875042516949852805536222144906

30) 95.0292318080446952680509981871741

31) 98.1709507307907819735377591608520

ш = 280 к = 1; собственные значения:

1) 3.83170597020751231561443588630041

2) 7.01558666981561875353704998148102

3) 10.1734681350627220771857117767752

4) 13.3236919363142230323936841269496

5) 16.4706300508776328125524604709880

6) 19.6158585104682420211250658841390

7) 22.7600843805927718980530051521805

8) 25.9036720876183826254958554459816

9) 29.0468285349168550666478198835296

10) 32.1896799109744036266229841044623

11) 35.3323075500838651026344790225187

12) 38.4747662347716151120521975577171

13) 41.6170942128144508858635168050603

14) 44.7593189976528217327793527132127

15) 47.9014608871854471212740087225063

16) 51.0435351835715094687330346332256

17) 54.1855536410613205320999662145325

18) 57.3275254379010107450905042437520

19) 60.4694578453474915593987498083827

20) 63.6113566984812326310397624178735

21) 66.7532267340984934153052597500420

22) 69.8950718374957739697305364354993

23) 73.0368952255738348265061175690917

24) 76.1786995846414575728526146235340

25) 79.3204871754762993911844848724880

26) 82.4622599143735564539866106487806

27) 85.6040194363502309659494254933803

28) 88.7457671449263069037359164348535

29) 91.8875042516949852805536222144899

30) 95.0292318080446952680509981871736

31) 98.1709507307907819735377591608515

Таким образом, имеется совпадение с 30 знаками после запятой. Отметим, что с таблицей VI работы [5] расхождений нет, но приведённые результаты немного точней.

3. Пример применения функций Бесселя целого индекса

В произвольной области Г € Я2 с достаточно гладкой границей дГ рассмотрим задачу (7), (8):

Ди

и\дГ = °-

Здесь функция /(г) либо задана, либо /(г) = [д(г) + Ар(г)]и(г), где д(г) и р(г) — заданные функции, и в этом случае имеем задачу на собственные значения для оператора Лапласа. В дальнейшем будем считать, что /, д и р

— гладкие функции.

Пусть г = ф((), \С\ ^ 1 — конформное отображение единичного круга на область Г; тогда в плоскости ( формально получаем те же соотношения (7)-(8), где, однако, вместо и(г) и /(г) следует писать и(() = и(г(()) и

\Ф'«)\2/(г(С)).

Дискретизация нулевого уравнения Бесселя рассмотрена в [4]. Оказывается, что этого достаточно, чтобы построить дискретизацию двумерной

задачи для уравнения Лапласа [6]. В круге радиуса 1 дискретный оператор Лапласа имеет вид:

2 п

Н = N^21 Лк  Нк,

к=0

где штрих у знака суммы означает, что слагаемое при к = 0 берется с коэффициентом 1/2, Лк, к = 0,1,... ,п — матрица размера т х т (матрица к-го дискретного уравнение Бесселя); Нк, к = 0,1,... ,п — матрица размера N х N: Нк^ = ео8[к2п(г — ] )/N)], г] = 1, 2,... N, через  обозначено кро-некерово произведение матриц. Здесь по в выбирается равномерная сетка: вг = 2пl/N, N = 2п + 1, I = 0,1,..., 2п. По г сетка может быть произвольна,

в данном случае выбираем ти = (1 + еов(^ — 1)п/2т))/2, V = 1, 2,..., т. То есть в единичном круге выбираем т окружностей, и на каждой окружности равномерную сетку по в:

Ло = В-1,

где В — матрица дискретного оператора (обратного к оператору Бесселя), построенная в [6];

Лд = Л0 + 2 К, к = 1,2,...,п;

К — диагональная матрица с числами (1/г^ )2, V = 1, 2,... ,т на диагонали.

В произвольной области матрица дискретного оператора Лапласа имеет вид: Z-1H, где Z — диагональная матрица с числами zvl= \ф'(£^)|2, (иі = = ти ехр(ів1), V = 1, 2,... ,т; I = 0,1,..., 2п [6].

Примеры численных расчётов. Для аналитически заданного конформного отображения вычисления проводились для эпитрохоиды, т.е. области получающейся из круга конформным отображением г = ^ (1 + є$Пр), \^\ ^ 1, є < 1/(пр + 1) (в программе обозначаются ЕРБ1 и №). На сетке т=50, N = 61 получено 73 первых собственных значения (см. ниже) для є = 1/6, пр = 4. Жирным шрифтом выделены знаки, совпавшие с расчетными на сетке 30 х 41:

М = 50 N = 61; ЕРБ1 = 0.1(6); NP = 4; время счета около 12 часов;

собственные значения:

1) 2.38444650949604970517317500817529 27

2) 3.73481160264336270303800592942815 28

3) 3.73481160264336270303800592946108 29

4) 4.60299170636341650906486673149750 30

5) 5.21305408447636594415167361552140 31

6) 5.40987176983340399646015402762260 32

7) 5.95284020254957074266667664726953 33

8) 5.95284020254957074266667664735062 34

9) 6.85046245675845252009714841665951 35

10) 6.85046245675845252009714841671300 36

11) 7.13745065059079443266158516815009 37

12) 7.25964235762084739781928108470845 38

13) 7.43089698228522450380843999247177 39

14) 8.20641593737488872010002074487221 40

15) 8.46703033293640948409921693554772 41

16) 8.46703033293640948409921696650293 42

17) 8.64626933595839253986265920247716 43

18) 8.78022273361123217092732024702811 44

19) 8.78022273361123217092732073102321 45

20) 9.43239467155670143601071573004845 46

21) 9.88298891926391094289884123527280 47

22) 9.88298891926391094290885516854269 48

23) 9.96025412843765989976019019427798 49

24) 9.96067883497458249296432391828820 50

25) 10.0677169120140574151369945503258 51

26) 10.3810329681107550983274521904108 52

10.9363318197115286486061588596269

10.9363318197115286490566933688893

11.1612248721630717492558347679048

11.2703687840812923758623562592516

11.2703687840812924379441899497323

11.7170723532432755105961602064307

11.8116162760450746824244865820387

11.8116162760450750359238289782441

11.8373790709381441799296422710614

12.2479455414825741525876370390611

12.4578735789763446630870262555313

12.4840161602227866181362114410024

12.7174707310530584401839103924169

12.7174707310530591099976162945176

13.1733648658278880535709423991943

13.2265396248205802726012953547913

13.2576192483420038095986918397910

13.3763092234868143064268852369909

13.3763092234873412793268306810793

13.7284082241398722372390888657113

13.7284082241405295639937794836824

14.0845624387645828041014741776930

14.2896890432957022566899381892198

14.3468569599344682193714374693753

14.3660088868535119687103057007520

14.3660088870628754963926158794226

53 14.5260032444333450006408510849440 64 16.0388925568050805309490904588994

54 14.7823390227650328822139182427564 65 16.1932121816405151627558326565322

55 14.7823390227701605589790191475342 66 16.3265204999443634036312371287739

56 15.0607377174290784301671781394397 67 16.3265205609221678113712393690573

57 15.1087606415694495735270739799535 68 16.3429023570436304015647080529037

58 15.3618732658649948819304045893759 69 16.6567315100655668236466221185369

59 15.4412325849722392654167903907534 70 16.8040692371885309000559415659651

60 15.4412325852883705119617736204168 71 16.8338073184134952489995022549795

61 15.6236132210873097648998985968366 72 16.8914598929759515336473804892064

62 15.7224016836042896639936517952767 73 16.8914602236018460187330833459450

63 15.7224017024130703061339497646196

Обсуждение полученных результатов. Расчёты производились на ПЭВМ Pentium IV с тактовой частотой 3,00 ГГц и объемом оперативной памяти 1 Ггб. Время последнего расчёта на сетке 50 x б1 около 12 часов. Предыдущий расчёт на сетке 30 x 41 занимает около 30 минут. Как видно из сравнения результатов на двух последних сетках, надежно с б-7 знаками после запятой определяется 30 собственных значений. Первое собственное значение определяется с 21 знаком после запятой. Расчёты проводились с учетверенной точностью REAL*^ с использованием транслятора с Фортрана Intel Visual Fortran 9.1. В приведенной выше таблице оставлены все знаки, которые даёт учетверенная точность. Знакам, выделенным жирным шрифтом, можно определенно верить. Первые б собственных значений получены на сетке 30 x 41 на ЭВМ БЭСМ-б [б, таблица 3.2, стр. 48] . Все знаки, приведённые в этой таблице, совпали (кроме последней цифры, по которой проводилось округление).

4. Нули функций Бесселя полуцелого индекса Jn+1/2

Вначале приведём результаты вычислений в Maple 11 для n = 0,..., 10: eval^BesselJZeros^^,! ..10));

3.141Б92бБ4, б.2В31ВБ307, 9.4247779б1, 12.Ббб370б1, 1Б.7079б327, 1В.В49БББ92, 21.99114ВБВ, 2Б.13274123, 2В.274333ВВ, 31.41Б92бБ4

evalf (BesselJZeros^^, 1 ..10));

4.4934094БВ, 7.72Б2Б1В37, 10.904121бб, 14.0бб19391, 17.2207ББ27, 20.3713029б, 23.Б194Б2Б0, 2б.ббб0Б42б, 29.В11Б9В79, 32.9Бб3В904

evalf (BesselJZeros^^, 1 ..10));

Б.7б34Б9197, 9.09Б011330, 12.32294097, 1Б.Б14б0301, 1В.бВ903б3б, 21.ВБ3В7422, 2Б.012В0320, 2В.1б7В2971, 31.32014171, 34.47048833

evalf (BesselJZeros(3.Б, 1 ..10));

б.9В7932000, 10.41711ВББ, 13.б9В0231Б, 1б.923б2129, 20.121В0б17, 23.30424б99, 2б.47б7б3бб, 29.б42б04Б4, 32.В0373239, 3Б.9б140БВ0

evalf (BesselJZeros(4.Б, 1 ..10));

В.1В2Бб14Б3, 11.7049071Б, 1Б.039бб471, 1В.3012ББ9б, 21.Б2Б41773, 24.727БбБББ,

27.91ББ7б20, 31.09393321, 34.2бБ39009, 37.43173б77

eval^BesselJZeros^^, 1 ..10)); 9.3ББВ12111, 12.9ббБ3017, 1б.3Б4709б4

29.332Бб2БВ, 32.Б24бб129

evalf (BesselJZeros^^, 1 10.Б12В3Б41, 14.2073924б 30.7303В073, 33.93710В30

evalf (BesselJZeros^^, 1

11.бБ703219, 1Б.4312В921 32.11119б24, 3Б.3331941В

evalf (BesselJZeros^^,! . 12.79078171, 1б.б41002ВВ 33.47бВ00В2, 3б.714Б2913

evalf (BesselJZeros(9.Б, 1 13.91Б822б1, 17.838б4320 34.828б9бБ4, 38.08247909

evalf (BesselJZeros(1С.Б, 1 1Б.0334б930, 19.02Б8Б3Б4 3б.1б81Б714, 39.43821448

3Б.707Б7б9Б

.10));

17.б4797487

37.13233172

.10));

18.92299920

38.Б413б48Б

10));

20.1824707б

39.93б12781

.10));

21.42848б97

41.3178б4б9

..10));

22.бб2720бб

42.б87бБ128

Б34б3411, 32

Ниже приводятся расчёты, аналогичные пункту 3: MACHEP =

19.бБ31Б210,

38.883б309б

22.904ББ0бБ, 2б.1277Б014,

20.9834б307, 24.2б27б804, 27.Б078б83б, 40.318892Б1

22.29Б34802, 2Б.б028ББ9Б, 28.8703733Б, 41.7390Б287

23.Б9127482, 2б.92704078, 30.2172б271, 43.14Б42Б02

24.87321392, 28.2371343б, 31.ББ018838, 44.Б39144б3

2б.1427б7б4, 29. 4Б.9212017б

870Б34б0,

1.9E-34.

m = 320, п=0.Б, л/ХЇ, i = 1,..., 10:

1) 3.14159265358968148253314596484760

2) 6.28318530717980998847751005721168

3) 9.42477796076904444882659617826803

4) 12.5663706143596199745004153513290

5) 15.7079632679484074188026742442073

6) 18.8495559215394299556118042407448

7) 21.9911485751277703949203122405416

8) 25.1327412287192399293501443109994

9) 28.2743338823071333796442247984913

10) 31.4159265358990498932469544081448

m = б40, п=0.Б, л/ЛЇ, i = 1,..., 10:

1) 3.14159265358978974602450302670173

2) 6.28318530717959346180037664523438

3) 9.42477796076936923807827248853378

4) 12.5663706143591869235912261908801

5) 15.7079632679489487301463333079038

6) 18.8495559215387803853630192704717

7) 21.9911485751285282222382168356239

8) 25.1327412287183738471062219818291

9) 28.2743338823081077143634600934317

10) 31.4159265358979673088112936148603

Таким образом, получается, что квадратные корни из собственных значений кратны п. Для сравнения приводим значение п с 4Б знаками, посчитанное в Maple 11: 3.141Б92бБ3Б897932384б2б4338 32 79Б028841971б940. Следовательно, получить больше 1Б верных знаков после запятой не получается. Задача с полуцелым индексом более сложная для вычислений. Дальнейшие вычисления приведены ниже. Отметим, что для n = 1, 5 проводились расчёты с m = 2240, в первом нуле получено 23 знака после запятой, а в 31 — 20 знаков после запятой.

5. Пример применения функций Бесселя полуцелого индекса

Рассмотрим задачу о собственных колебаниях сферы радиуса г0 с нулевыми граничными условиями первого рода. Эта задача сводится к отысканию собственных значений и собственных функций уравнения Аь + Хь = = 0 с граничным условием на поверхности сферы V = 0, [7]. Обозначим

(га) (га) (га) т ґ \ г\

VI , V2 ,...^т корни трансцендентного уравнения ^га+1/2^) = 0, нахо-

дим собственные значения Хт,п = . Каждому собственному значе-

нию Хп,т соответствует 2п+1 собственная функция. Введём обозначение

Вычисления проводились с помощью подпрограмм: QBESSEL вариант с учетверенной точностью программы BESSEL [4,6], QTELM фортранный вариант с учетверенной точностью алгольной программ elmhes [8], N — размер матрицы, H — исходная матрица заменяется на выходе матрицей в верхней форме Хессенберга; QHQR фортранный вариант с учетверенной точностью алгольной программ hqr [8], N — размер матрицы, H — исходная матрица в верхней форме Хессенберга (не сохраняется), MACHEP — машинная е, равная наименьшему положительному числу для которого 1 + е > 1, WR, WI выходные массивы длины N, которые содержат на выходе действительную и мнимую часть вычисленных собственных значений, CNT — целый массив длины N, который содержит на выходе число итераций для вычисления каждого собственного значения. QMINV — вариант с учетверенной точностью подпрограммы MINV для обращения матрицы [9].

Заключение. По поводу получения полных версий программ обращайтесь по электронному адресу: algazinsd@mail.ru или на адрес Института проблем механики РАН, 119526, Москва, проспект Вернадского д.101, к.1.

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

2. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 744 с.

фп(х) = \[П^п+1 /2(х). Тогда собственные функции, рассматриваемой задачи, можно представить в виде:

(п = 0,1,m = 1, 2,j = -n,..., -1, 0,1,..., n), уП0)(в,ф) = Pn(cos в), уП-^(е,ф) = pnj)(cos в) cos jф,

Уп](в,Ф) = Pij)(cos в) sin jф, (j = 1,2, ...,n).

Pnj — присоединенная функция Лежандра j-го порядка [7]. Примечание.

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

3. Гончаров В.Л. Теория интерполирования и приближения функций. М.: Гостех-теориздат, 1954. 328 с.

4. Алгазин С.Д., Бабенко К.И., Косоруков А.Л. О численном решении задачи на собственные значения. Препринт. М.: ИПМ, 1975. №108. 57 с.

5. Таблицы нулей функции Бесселя. Библиотека математических таблиц. М.: ВЦ АН СССР, 1967. Вып.44. 95 с.

6. Алгазин С.Д. Численные алгоритмы классической математической физики. М.: Диалог-МИФИ, 2010. 240 с.

7. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 742 с.

8. Уилкинсон Дж, Райнш Ц. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. М.: Машиностроение, 1976. 391 с.

9. Сборник научных программ на фортране. Матричная алгебра и линейная алгебра. М.: Статистика, 1974. Вып.2.

Алгазин Сергей Дмитриевич (algazinsd@mail.ru), д.ф.-м.н., ведущий научный сотрудник, Институт проблем механики им. А.Ю. Ишлинского РАН, Москва.

About tabulation with split-hair accuracy of zero of cylindrical

functions

S. D. Algazin

Abstract. The algorithm of calculation of zero of cylindrical functions which allows calculating some first tens zero from 15-30 signs after a comma is reduced.

Keywords: Bessel equation, numerical algorithms without saturation.

Algazin Sergey (algazinsd@mail.ru), doctor of physical and mathematical sciences, leading researcher, Institute for Problems in Mechanics of RAS, Moscow.

Поступила 17.01.2013