Алгоритм 20а получен в результате исправления, сокращения и ординарной переработки алгоритма 20 (Peavy S. «САСМ», 1960, № 10). Кроме исправления опечатки, указанной в соответствующем «Замечании» CPeavy S. «САСМ», 1961, № 2) и «Подтверждении» (Alexander W. J., Thacher Н. С «САСМ», 1961, №4), в алгоритме 20 была исправлена опечатка- отсутствие точки с запятой перед последним оператором тела процедуры.

Подтверждение к алгоритму 20

В. Дж. Александер и Г. Тачер

(Alexander W. J., Thacher Н. С. «САСМ», 1961,

№ 4)

Процедура ехрШ{х) была запрограммирована для машины LGP-30, использующей 7-разрядный т1ранслятор с плавающей запятой и 8-раэрядную интерпретирующую программу 24.2. Константы, заданные более чем семью значащими цифрами (нли более чем восьмью цифрами для программы 24.2), соответственно округлялись.

После исправления константы 0.005519968 на 0.05519968 обе црограммы дали удовлетворительную точность В проверяемой области.

Результаты, полученные по 8-разрядной программе 24.2, сравнивались с данными 9-разрядных математических таблиц [211] для последовательности значений х=0.1, 1.0 и 10.0. Наибольшее расхождение обнаружено для x=0.1 и равно -16X110-*. Для х>1 все проверяемые значения оказались точными до восьми разрядов включительно.

При вычислении вещественных значений показательной интегральной функции этот алгоритм работает быстрее, чем ekz (алгоритм 14). Для x<i\ отношение скоростей было порядка 20.

Свидетельство к алгоритму 216

Алгоритм 216 не публикуется здесь потому, что позднее в журнале «САСМ» для вычисления бесселевых функций первого рода (и не только для щелых порядков)

был .опубликован более совершенный алгоритм 236 («САСМ», 1964, № 8). С;м. также алгоритм 236а [28].

АЛГОРИТМ 226

Функция Риккати-Бесселя первого и второго рода

[S17]

Эта процедура вычисляет функции

для вещественного хф{) и всех целых значений ft от О до п с заданной (абсолютной) точностью eps. Вычисление производится с использованием рекуррентных соотношений для цилиндрических функций. Для n<abs{x) как Sfe(x), так и Си{х) вычисляются с использованием рекурсий в возрастающих порядках. Для n>abs{x) функции Sh{x) получаются с использованием рекурсии в убывающих порядках [191, с. 255-257].

Два вектора S (х) и S (х), взятые из двух различных интервалов порядка выше п, проверяются, удовлетворяет ли максимальная компонента их разности допустимой погрешности eps. Если не удовлетворяет, то для достижения требуемой точност.и проводится максимум 10 итераций. Начальные значения Skmax и Shmax-i ДЛЯ обратной итерации определяются по соответствующим значениям Chmax~i И Chmax- В случэе n<abs{x) никакой проверки точности «е делается. В этом случае как Ck{x), так и Sh(x) могут получить ошибки того же порядка, что и подпрограммы для sin {х) и cos (х).

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

При х==0 - выход к глобальной метке signal22. procedure riccatibessel (x,n,eps,acc,stp,imax) result: (s,c); value x,n,eps,acc,stp,imax; real x,eps,acc,stp,imax; Integer n; array s,c;

Ucmn real rl,r2,r3,r4,r5,r6,a,b,dl,cl2,p,max; integer i,k,m; array MO:n];

procedure гг(с,р,г1,г2,гЗ); value p; real р,г1,г2,гЗ; array c; begin c[p]:=rl:=(2xk-l)Xr3/x-r2;

г2: = гЗ; r3: = rl end rr;

start: i: = l;

if x=0 then go to signal22; if n<abs(x) then easel: begin rl:=-sin(x);

r2:=г4: = cIO]: = cos (x);

r5: = sI0]: = sin(x);

for k: == 1 step 1 until n do

begin rr(c,k,r3,rl,r2); rr(s,k,r6,r4,r5) end; go to fin end easel; case2; m: = l; rl:=-sin(x); r2:=ci0]: = cos(x);

for k:= 1 step 1 until n do гг(с,к,гЗ,г1,г2); a:=n;

iter: p:= if abs(xXll then 12 + a else 2Xa + l; for k:= n + 1 step 1 until p do begin dl:=(2Xk-l)/x; if ln(abs(dl))+ln(abs(r2))>18 then go to signal22; r3: = dlXr2-rl;

if abs(r3Mn])>acc then go to out; rl: = r2; г2: = гЗ;

comment Этот цикл наиболее подвержен

переполнениям

end iter; г2:=г1; к:=р;

г6:=0.25Х(х/к)2/г2; г5: = 1/гЗ; if m=2 then go to norm; m:=2;

for k:= к step -1 until 2 do

rr(w, if k>n + 2 then n else k~2, г4,г6,г5); dl:=-r5/x-г6;

d2:= if abs(w[0])>abs(dl) then sin(x)/w[0] else cos(x)/dl;

for k:= 0 step 1 until n do w[k]: = d2Xw[k]; acc:=stpXacc; a: = a+stpf (1/3); r2: = c[n]; rl:=ci;n-1];,