Главная страница Алгоритмы [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
Эти расчеты проводились для 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.013 |