К разделу 'Математика' сайта AlgoList.

 
Это первый документ - нажмите "Оглавление"ОглавлениеВперед

Основные разложения

В этом параграфе описаны алгоритмы, позволяющие, в принципе, вычислять norm_x с произвольной точностью: ошибки всегда возникают из-за конечной разрядности. Алгоритмы основаны на разложении в ряды Тейлора и цепные дроби функций, связанных с norm_x простыми соотношениями. Структура всех этих алгоритмов описывается в приложении А; коды на Си содержатся в приложениях В и Г.

Вероятно, наиболее известными соотношениями, позволяющими вычислять norm_x с произвольной точностью, являются
 

norm_x = 0.5 + I(x),
(1)
.

Это разложение легко получить, интегрируя почленно разложение Тейлора в нуле функции exp(-x2/2). Хотя при этом получается знакопеременный ряд, первый отброшенный его член не является, вообще говоря, верхней границей для остатка, т.к. его члены монотонно убывают, лишь начиная с некоторого. Тем не менее, довольно легко убедиться в том, что при NrmDF03=0.3989... это приятное свойство имеет место. А так как в любой практической задаче eps <<0.1, можно считать, что если выполнено неравенство

(2) NrmDF06,

то сумма n первых членов ряда доставляет I(x) с точностью eps (если, конечно, не учитывать накопление ошибки).

Алгоритм I. Этот алгоритм основывается на разложении (1). Обратите внимание: при n-м выполнении шага I2 переменная m равна 2n, так что t здесь делится все-таки на 2nn!. Разложение (1) было использовано в алгоритме 272 [12], который, однако, значительно менее экономен, чем нижеследующая версия.

I1.Положить t=sum=x; x2=x*x; m = 2.
I2.Положить t=-x2*t/m; s=sum; sum=s+t/(m+1); m=m+2.
I3.Если |s-sum| > e, то перейти к I2.
I4.Положить F = 0.5 + sum/Ц2p.
Менее известными, но столь же просто получаемыми, как (1), являются соотношения (3) , где .

Это разложение для можно получить почленным интегрированием ряда Тейлора функции j(x)/j(t) в нуле. Для остатка ряда (3) справедливы оценки

, поэтому оценка Fnошибки ограничения для F(x) удовлетворяет неравенству

(4) .

Неравенства (2) и (4) показывают, что скорость сходимости разложений (1) и (3) падает с ростом x.

Алгоритм T. Этот алгоритм основан на разложении (3). В работах [1, 11] в аналогичном алгоритме на шаге T1 производится присвоение t = sum = |x|*exp(-x2/2)/Ц2pс очевидным изменением шага T4. Так как, однако, шаг T3 там не изменен, условием останова оказывается |sn+1 - sn| < e*exp(x2/2)/Ц2p, что, согласно (4), не гарантирует необходимую точность вычисления F(x). Тем не менее, если вычисление F(x)производится с машинной точностью (и машинный нуль достаточно мал) и ошибка при вычислении экспоненты не слишком велика, алгоритмы в [1, 11] дают в пределах представления чисел с плавающей запятой правильный результат.

T1. Положить t = sum = |x|; x2 = x*x; n = 3.
T2. Положить t = x2*t/n; s = sum; sum = s+t; n = n+2.
T3. Если |sum - s| > e , перейти к T2.
T4. Положить F = 0.5 + sign(x)*sum*exp(-x2/2)/Ц2p .
Вот еще одно разложение функции , которое можно применять для вычисления F(x)при небольших x:

(5) 

Эту цепную дробь можно получить разными способами (см., например, [2,3]). Члены ее звеньев альтернируют и потому характер сходимости ее подходящих дробей довольно сложен. Однако, для подпоследовательности подходящих дробей с четными номерами можно показать, что ее четные члены сходятся к монотонно возрастая, а нечетные – монотонно убывая. Поэтому модуль разности |w2n-2-w2n| двух соседних подходящих дробей этой подпоследовательности оценивает сверху уклонение каждой из них от .

Алгоритм S. В данном алгоритме, основанном на разложении (5), для получения подпоследовательности подходящих дробей с четными номерами мы производим на каждом витке цикла вычисление двух очередных подходящих дробей; соответствующее рекурсивное соотношение для сжатой дроби слишком сложно и потому не использовано.

S1. Положить x2 = x*x;
S2. /* Вычисление подходящих дробей с номерами 0 и 2:
        предпервое состояние цикла */
     rho = 1/(3-x2); sum = |x|*rho;
     term = x2*sum; rho=3*rho; s = sum = 3*sum;
     n = 1;
S3. /* Вычисление четной (когда x2 < 0)
        или нечетной (x2 > 0) подходящей дроби */
     Положить n = n+1; r = x2*m/(4*n*n-1);
     rho=1/(1+r*rho); term=(rho-1)*term;
     t = s; s=sum; sum=sum+term;
     x2=-x2;
S4. Если x2 < 0 (т.е. на шаге S3вычислялась нечетная
    подходящая дробь), то перейти к шагу S3.
S5. /* Проверка точности */
   Если |s-t| > e или |s-sum| > e, перейти к S3.
S6. Положить F = 0.5 + sign(x)*sum*exp(-x2/2)/Ц2p;
     на этом работа алгоритма заканчивается.
Соотношения [13]

(6') wn-1-wn = (-1)n(4n-1)(2n-2)! x4n-1/(q2n-2q2n),

(6") 

показывают, что необходимое число подходящих дробей растет с ростом x. Поэтому разложение следует применять лишь при небольших значениях x; сам Л.Шентон [13] рекомендует применять его только при x <Ц3.

До сих пор мы имели дело с разложениями, "качество" которых ухудшалось с ростом x. Рассмотрим теперь соотношения, которые с ростом x позволяют вычислять F(x) все быстрее:

(7') F(x) = 0.5 + j(x)R(x),

(7") 

Функция R(x) называется отношением Миллса. Разложение в цепную дробь для него было предложено еще Лапласом в 1805 г. Члены звеньев этой цепной дроби при x > 0 положительны и потому (см., например, [6]) подходящие дроби с четными номерами сходятся к R(x) монотонно возрастая, а с нечетными – монотонно убывая. Поэтому модуль разности

двух соседних подходящих дробей этой подпоследовательности оценивает сверху уклонение каждой из них от R(x). Можно показать, что имеет место оценка

(8) ,

из которой следует, что скорость сходимости разложения растет с ростом x.

Алгоритм L. В этом алгоритме F(x)вычисляется с помощью разложения отношения Миллса R(x) в цепную дробь Лапласа (7).

L1. Положить x2=x*x; y=exp(-x2/2)/Ц2 p; x2=1/x2;
      sum=term=y/|x|; r=s=0; rho=1.
L2. /* Вычисление очередной подходящей дроби */
    Положить r=r+x2; rho=1/(1+r*rho); term=(rho-1)*term;
      t=s; s=sum; sum=sum+term.
L3. /* Проверка точности */
    Если |s-t| > e или |s-sum| > e , перейти к L2.
L4. Если x > 0, положить F = 1-sum и остановить алгоритм.
L5. Положить F = sum.
Для вычисления F(x) при больших x используют также асимптотический ряд

(9) ,

который легко получить, проинтегрировав отношение Миллса R(x) по частям. По заданному e можно найти xe, такой, что при x і xe, ряд (9) позволяет вычислять F(x) с точностью e. Однако, с помощью этого ряда невозможно, вообще говоря, достичь в любой заданной точке x произвольной точности, поэтому мы и не исследуем алгоритм, основанный на этом разложении.
 
 

Для проведения численных экспериментов алгоритмы I, T, S и L были реализованы на Си. В нижеследующей таблице для каждого алгоритма в столбце N приводится количество "витков" цикла, потребовавшихся для вычисления значения с e = 0 для нескольких x; поскольку все алгоритмы дали одинаковые 8 знаков после запятой, сами вычисленные значения опущены. В столбце O приводится число потребовавшихся операций. Для алгоритма I оно вычислялось по формуле 8+11ґN, для алгоритма T – по формуле 14+8ґN, для S – по формуле 23+25ґN, наконец, для L – по формуле 16+10ґN. Обратите внимание – всем операциям, в том числе, таким, как вызов функции fabs или exp, придается равный вес. Для иллюстративных целей такая точность признана достаточной.
 
x
I
T
S
L
 
N
O
N
O
N
O
N
O
0.25
8
96
8
78
8
223
6027
60286
0.50
10
118
10
94
10
273
1507
15086
0.75
12
140
13
118
12
323
694
6956
1.00
14
162
15
134
14
373
391
3926
1.25
16
184
17
150
14
373
259
2606
1.50
18
206
19
166
16
423
180
1816
1.75
20
228
21
182
18
473
137
1386
2.00
23
261
23
198
20
523
108
1096
2.25
25
283
25
214
22
573
88
896
2.50
27
305
27
230
22
573
74
756
2.75
30
338
29
246
24
623
63
646
3.00
33
371
32
270
26
673
54
556
3.25
35
393
34
286
28
723
48
496
3.50
38
426
36
302
30
773
44
456
3.75
41
459
39
326
32
823
39
406
4.00
44
492
41
342
34
873
35
366
4.25
48
536
44
366
36
923
33
346
4.50
51
569
47
390
38
973
31
326
4.75
55
613
49
406
40
1023
28
296
5.00
58
646
52
430
42
1073
26
276
5.25
62
690
54
446
44
1123
25
266
5.50
66
734
57
470
48
1223
24
256
5.75
70
778
60
494
50
1273
22
236
6.00
74
822
63
518
52
1323
21
226
6.25
79
877
66
542
54
1373
21
226
6.50
83
921
69
566
58
1473
20
216
6.75
88
976
72
590
60
1523
19
206
7.00
93
1031
76
622
64
1623
18
196
7.25
98
1086
79
646
66
1673
18
196
7.50
103
1141
82
670
70
1773
17
186
7.75
108
1196
86
702
72
1823
17
186

Построим графики, отражающие данные из этой таблицы, – их проще анализировать. На нижеследующих рисунках слева – столбцы N, справа – O.
 

а) число витков цикла

б) количество операций

Мы видим, что по числу N витков цикла алгоритмы при не слишком больших x практически не различаются до x»3.5, после чего лидирует алгоритм S. Начиная с x» 4, бесспорным лидером становится алгоритм L. По числу операций, конечно, более важной характеристике, до x» 4 лидером является алгоритм T, после – снова алгоритм L.

Помимо превосходства в скорости алгоритм T выгодно отличается от I и S устойчивостью по отношению к накоплению ошибки. Для доказательства высказанного утверждения оценим распространяемую ошибку округления.

Пусть – относительная ошибка n-й частичной суммы в алгоритме I, а – относительная ошибка tn. Для простоты предполагаем, что представление целых чисел и действия над ними осуществляются без ошибок, а ошибки округления всех арифметических операций над действительными числами не превосходят r. Тогда для имеем следующее соотношение

(10) ,

причем

(11) 

где – относительная ошибка представления x2. Будем считать также, что

. Из (11) с легкостью следует , а из (10) получаем , откуда

(12) .

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

Для алгоритма T относительная ошибка частичной суммы sn также удовлетворяет соотношению (10), однако теперь . Действуя, как раньше, получаем оценку

(13) .

Эта оценка возрастает лишь немного медленнее, чем оценка (12), Однако, члены ряда (6) положительны, поэтому истинная относительная ошибка значительно меньше полученной оценки (13). Хотя общий член ряда (2) убывает быстрее, чем общий член ряда (3), машинные эксперименты показывают, что число исполнений шага T2 меньше числа «витков цикла» в алгоритме I; это также приводит к меньшей распространяемой ошибке. Тем не менее, на точность вычисления F(x) сильно влияет точность вычисления exp(-x2/2): возникающая здесь ошибка может «съесть» все преимущества разложения (3).

Полученные данные позволяют рекомендовать для вычисления F(x) с произвольной точностью сочетание алгоритмов T и L. Нужно только выбрать точку p, разделяющую области действия этих алгоритмов.

Итак, при возрастающих x время счета по алгоритму T возрастает, а по алгоритму L убывает. Поэтому в качестве p следовало бы выбирать точку, в которой эти времена становятся равными. К сожалению такая точка зависит от точности счета e , от разрядности используемого компьютера, от его скорости. Поэтому выбрать навеки оптимальную» точку p не удастся; я рекомендую выбирать ее при e , равном машинному нулю.
Вверх
Это первый документ - нажмите "Оглавление"ОглавлениеВперед

На правах рекламы
По хорошей цене наркологический центр в москве на любых условиях.