Главная страница  Алгоритмы 

[0] [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [ 15 ] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

go to iter; - r- .....

norm: for k:=: к step-1 until 2 do

rr(s, if k>n+2 then n else k-2,r4,r6,r5); dl: = r5/x-r6;

d2:= if abs(sIO])abs(dl) then sin(x)/s[0] else

cos(x)/dl; max:=0;

for k:= 0 step 1 until n do begin sIk]: = d2Xs[k]; b: = abs(sM-wM);

if b>max then max: = b end;

if max<eps then go to fin; for k:= 0 step 1 until n do wi[k]: = slk]; acc: = stpXacc; if imax then go to fin; i: = i+l; a: = a + stpt(l/3); r2:=cln]; rl: = cln-Ij; go to iter; tin: end riccatibessel;

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

Алгоритм 226 получен ;из алгоритма 22а путем внесения в него поправки, указанной Б. Томасом в нижеследующем «Подтверждении к алгоритму 22», а также внесения двух следующих модификаций.

1. Третья строка сверху на странице 36 [24]

begin гЗ:- (2хк-1) Хг2/х-rl; была заменена на

begin dl:=(2Xk-1)/х; if ln(abs(dl))+ln(abs(r2))>18 then go to

signal22; r3: = dlXr2-rl; для более точного определения момента и причины возникновения переполнения.

2. Оператор

out:r6:x \ 2/(4ха! f 2Х-2);

был заменен на оператор

ои1:гЬ: = 0.2БХИк) 2/г2;

для упрощения процедуры и экономии машинного времени.



Алгоритм 226 был транслирован «а машине М-220 в системе ТА-1М, и для него были повторены расчеты, проведенные ранее для -процедуры пеитап алгоритма 22а [24, с. 37-38]. Результаты сведены в табл. 9.

Таблица 9

Функция

5.(х)

влх) SAx)

0.14112Ю0

1.0370324

0.89591249

-0.27941549

-1.0067395

-0.22395427

0.41211848 0.95692120 -0.093144750

-0.53657291 -0.88856836 0.31443082

г. (X) с, (X) Ci (X)

-0.98999249 -0.18887749 0.80111500

0.96017028 -0.11938711 -1.0198638

-0.91113026 0.31088178 1.0147575

0.84385395 -0.46625175 -0.96041689

Эти расчеты проводились для acc=W, stp = l(fi, imax=lOO, п=2 и epslOr-\ 10- причем результаты для eps=10~* не отличались от"результатов для eps= =10-. Расхождение в знаках и значениях 52(9) из вышеприведенной табл. 9 по сравнению со значениями из табл. 9, приведенной в «Свидетельстве ,к алгоритму 22а» [24, с. 38], объясняется ошибками в последней, поскольку правильность вышеприведенной табл. 9 была подтверждена также н трансляцией процедуры пеитап (алгоритм 22а), результаты .которой полностью совпали с данными этой новой табл. 9..

Кроме того, были проведены расчеты для х=ОЛ, 0.5, 1, 2, 3, 5, 7.5, 10, 50 с целью определения максимального значения п, позволяющего вычислять с помощью процедуры riccatlbessel без переполнения. В результате этого расчета оказалось, что процедура riccatibessel позволяет вычислять без переполнения (по .крайней мере на машинах типа М-20) лишь при nentler{x)+ \ для нецелых X и при пх для целых х.

Для процедуры пеитап были проведены расчеты при д;=0.1, 5, 10 и п=5, 10, 15 соответственно. Эти расчеты показали, что процедура пеитап позволяет вычислять S и С без переполнения и при п>х. Результаты расчетов с помощью процедуры пеитап полностью (во всех выданных на печать восьми цифрах) совпали с результата-Ми расчетов с помощью процедуры riccatibessel. Таким образом, проведенные расчеты показали превосходство



йроцёдуры neuman над процедурой riccatibessel по размерам области применимости. Кроме того, превосходство процедуры пеитап по быстродействию очевидно.

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

Алгоритм 22а получен в результате ординарной переработки, некоторого сокращения и модификации алгоритма 2-2 (Oser Н. «САСМ», 1964, № 11).

В соответствии с конечными фррмулами, приведенными в [3, с. 202-203; 4, с. 26], было составлено описание процедуры пеитап, вычисляющей те же функции Риккати - Бесселя. При этом учитывались соотношения

(имеется виду, что lft+iy2(-) = (-1)~- й 1/2 W) Д-"» хфО \ik = 0, \,..:,п) [3, с. 203].

procedure neuman (x,n) result: (s,c);

value x,n; real x; integer n; array s,c; begin real y,y2,a,b,aO,bO,p,slI,sl2,s21,s22; integer k,j,kl; s[0]:=sin (x); d[0]:=cos (x); y: = 2Xx; у2:=уХу; for k:=.l step 1 until n do begin a: = aO: = l; b:=kx (k+l)/y; kl:=k-f-(2; for j:= 1 step 1 until kl do begin p:=2xj; sl 1:= (k-p + 1) X (k-p+2) X (k + p-1) X

(к+р)/(рХ(р-1)Ху2); s21:=(k + p + l)x(k-p)/(yX(p + l)); aO:=-sllXaO; a: = a + aO; b0:=:a0xs21; b: = b + bO end I;

j:=if kl = (kl2)x2 then 1 else -1; sll:=sl2: = s21:=s22:=jxs[0]; if к=И=к1х2 then sll: = s22:=-jXc[0] else begin sl2: = s21:==jXc[0]; b:=b-bO end; s[k]: = axsll+bxsl2; c[k]: = aXs21-bxs22 end к end neuman;

Процедура neuman была проверена с помощью транслятора ТА-1 для п=4 и х-Ъ, 6, 9, 12. Кроме того, были вручную запрограммированы формулы из [4 с. 26]




[0] [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [ 15 ] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

0.0299