Доставка цветов в Севастополе: SevCvety.ru
Главная -> Справочник по алгоритмам

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 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

вычисляется разложением в ряд

а {X)

f + \nx+J

2/г-(2«)!

Программа 6.6.

05 PRINTВЫЧИСЛЕНИЕ ФУНКЦИИ CliXJ

10 INPUT ВВЕДИТЕ Х=Х

20 LETS=X"£:LETI=0:LETA=0

30 LETB.=l:LETC=l

40 LET1=1+1:LETD=2*I

50 LETC=<D-1>*D*i:-C): LETB=B*S

60 LETE=B-- < СжП?: LETA=fi+E

70 IF ftBS<EKlE-9 THEN 90

80 БОТО 40

90 PRINT CI<X>=A+.5772156649+L0&«:X-> 100 C-.OTO 10: END

Пример. Для x= 1 получим Ci (1) = =0,337403923.

Вычисление Si (x) и Ci (x) численным интегрированием целесообразно лишь при малых значениях х. При больших х> 8 существенное уменьшение времени вычислений достигается при использовании асимптотических выражений [36].

§ 6.4. Гамма-функции (включая неполные)

Гамма-функция

Г (х)=Г (X) =5 е-/- dt

является одной из наиболее распространенных специальных функций. Поэтому целесообразно рассмотреть различные алгоритмы ее вычислений.

Для - 18<х<49 с погрешностью порядка 1 • 10-* Г (х) может быть вычислено с помощью 20-кратного преобразования по формуле Г (х+1)=хГ (х) и с применением

Программа 6.7.

затем формулы Стирлинга Г(х)=.д/

(6.1)

Такое преобразование необходимо для уменьшения погрешности (6.1). Алгоритм вычисления Г (х) при этом следующий.

1. Вводим X и полагаем Z=x, X = Z и А = \.

2. Вычисляем Z,-!-Z, (X + I), где / = = 1, 2, 20.

3. Полагаем В = Х--21 и вычисляем

,G = exp (B(lnB-l) + l/(12B))-\/Z2o. Программа 6.7.

Пример. Для различных х получаем Г ( - 1,5) =2,363272685 (2,363271801),

Г (- 0,5) = - 3,544908843 (- 3,544907702). Г (1,5) =0,8862271413 (0,886226926), Г (40) = =2,039788231 -10" (2,039788208- Ш**). В скобках указаны точные значения Г (х).

Более высокую точность при - 10<х<70 имеет метод, в основе которого лежит следующий алгоритм.

1. Вводим х=Д и задаем А=А - 1.

2. Полагаем Z=y4 -64.

3. Если Z-<0, идем к п. 4, иначе задаем £) = 1 и идем к п. 6.

4. Вычисляем B=lint (Z-1) и полагаем С=А, А=А + В я D=\.

5. Находим

С, = С, -Ь1, D, = D, ,/C,

где / = В, В-1, В -2, 2, 1.

6. Вычисляем G = Г (х) по формуле Стирлинга в следующей ее модификации:

G = (A) eO/(i2A))D

При этом алгоритме формула Стирлинга используется при аргументе, значения которого превышают 64, что дает меньшие погрешности. Однако при х<; -10 корректирующий множитель ZX1-10-** и ему, как и значению Г (х), приписывается нулевое значение.

Программа 6.8.

Пример. Г ( - 1,5) =2,363271825,

Г (- 0,5) = - 3,544907737, Г (1,5) =

=0,8862269342 н Г (40) =2,039788229-10«.

10 PRINTВЫЧИСЛЕНИЕ ГЙММЙ-ФУНКЦИИ Г<Х>

20 INPUTBBEAHTE X=2sLETX=2!LETft=l

30 FOR 1=1 ТО 20

40 LET2=Z»;«:x+I>sNEXT I

50 LETB=X+21

60 LETG=EXP(B»;(L06<B>-1 >+l.12.B>жЗ£!Р«:2ж#Р1В>2 70 PRINT.ЗНАЧЕНИЕ r<X)=6sG0T0 20: END

Программа 6.8.

10 PRINTВЫЧИСЛЕНИЕ ГАММА-ФУНКЦИИ ГСХ?

£0 INPUTВВЕДИТЕ X=X:LETA=X-l:LETZ=A-64:IF 2<0 THEH 48

30 LETD=l:&OTO 60

48 LETB=ABSaNT<Z-l)):LETC=A:LETA=A+B:LETD=l:FOR I=B TO 1 STEP -1 50 LETC=C+1:LETD=D.C:NEXT"I

60 LETC-.= i < A/-EXP < 1 > ) " A > ЖЕХР a 1 A > *SQR < 2*-#P I ж A > *D 78 PRINTr<X)=&:GOTO 28:END

a(x)=CI (X)=v + lnA:+t dt



Наименьшее время вычислений обеспечивает алгоритм, основанный на аппроксимации Г {х) на отрезке [О, 1] степенным полиномом.

1. Из аргумента B=x последовательно вычитаем 1 для того, чтобы В попал в отрезок [О, 1]. Одновременно подсчитываем произведение £) = В(В-1) (В -2) ... (В - - К) у где /( - число вычитаний 1 из В.

2. Если лг> О, вычисляем Г (х) по формулам Y=B~K. Fr(Y+\) = = (((((((0,035868343} - 0,193527818) Y + + (0,482199394) Y - 0,756704078) Y + -f 0,918206857) Y - 0,897056937) Y + + 0,988205891) Y - 0.577191652) К + 1, Г (X)=DF/X.

3. Если x<;0, вычисляем К и f по п. 2 для х1, после чего находим

sinnx Г(х-1) sinnx DF

Программа 6.9.

Пример. (1,5) ==0.03648992674.

11) (50) =3,901989672.

Неполные гамма-функции задаются соотношениями

Р{а,х)=-\е-Ч-Ш,

у (а. X) =Р (а, X) Г (а) = e-t"- dt.

(6.2)

Г (а, х) = Г (а) -y (а. х) =\ e-t"- dt, о

у* (а, х) =х°Р (а, x) =x-°y (а- (а).

В качестве основной можно выбрать функцию у{а,х), легко вычисляемую разложе-

10 PRINTВЫЧИСЛЕНИЕ ГЙММЙ-ФУНКЦИИ ГСХ> ПО АППРОКСИМАЦИИ

£0 1НРиТВВЕДИТЕ K=X:LETB=ABS«)!LETD=1

36 IF В<=1 БОТО 50

40 LETD=D*B:LETB=B-l!60T0 30

50 LETF=«:«:.035868343*B-. 1935£7818>жВ+.4821993Э4)жВ

66 LETF=< i<F-.7567846785*8+.Э18286857)жВ-.897856937>*B

70 LETF=< С F+.Э88£0589 О *B~,57719165£)*B+1

80 LET&=F*ri/X! IF y.<B GOTO 100

90 PRINTr«j=G:GOTO £0

lee LET&=#PI/SIN<#Pl!«X>/D/F!GOTO 90sEND

Пример. r (- 3,2) = 0.6890562856,

Г (- 2,5) = ~ 0,9453086453, Г (2) = 1,

Г (0,5) = 1,772453992 и Г (50) =6,082818641 X X 10 . Погрешность вычислений в этом методе не превышает З-Ю"".

Логарифмическая производная гамма-функции (диагамма-функция)

(х)=-[1пГ(х)]

с погрешностью менее 2-10~* для х> 20 может вычисляться по асимптотическому разложению тф (х) = I п х - 1 / (2х) -

- 1/(12х)+ ... С помощью рекуррентной формулы ij) (х -f 1) = ij) (х) + 1 /х значения ij) (х) при х<20 можно найти, вычисляя ij) (х + 20) и последовательно применяя эту формулу в направлении уменьшения х. Программа 6.10.

нием в ряд (х>0, а> 0) y (а, X) =

aLL. (a+\)(a + 2)...{a+i) + ew-Hl].

где yV = int (2х + 8) и B/. - корректирующий член. Зная у(а,х) и Г (а), можно вычислить любую другую из приведенных неполных гамма-функций. Вычисление v (а, х) реализуется программой 6.11. Программа 6.1 i.

Пример. y (2; 8) =0,9969825948,

y (0,5; 1) = 1,493648313.

Неполная гамма-функция у (а, х) может вычисляться также численным интегрирова-

10 PRINTВЫЧИСЛЕНИЕ ЛОГАРИФМИЧЕСКОЙ ПРОИЗВОДНОЙ

15 PRINTГАММА-ФУНКЦИИ IKLOG Г<Х))>11Х

£0 INPUTВВЕДИТЕ X=X5IF X>26 THEN 40

30 LETC=20+«:X-INT«:X>>sGOTO 50

40 LETX=X+l!LETC=X

50 LETA=l/2/-C! LETA=-A*A/3-A+L0G«:C)

60 LETC=C-l:LETA=A-l/C

70 IF X<C THEN 60

88 PR I NTЗНАЧЕНИЕ IKLOG Г<ХУ) Г- = (ч 90 GOTO 20:END

Программа 6.11.

16 PRINTВЫЧИСЛЕНИЕ НЕПОЛНОЙ ГАММА-ФУНКЦИИ G«:A,X> 20 INPUTВВЕДИТЕ A/X A,XsLETB=INTC2*X+9>sLETV=B 30 FOR I=B TO 1 STEP -lsLETV=V*Xa+A)+lsNEXT I 40 LETG=V*EXP«:-X>*X"A/A:PRINT&«:A/X>=G!G0T0 2@!END



нйем - см. формулу (6.2). Так, используя программу 4.42 численного интегрирования методом Гаусса и задав а = 2, В=х=8, Л=0 и M = 8, получим Y (а, х) =0.9969832028 при с»20 с.

§ 6.5. Функции Бесселя (включая модифицированные)

функции Бесселя являются решениями линейного дифференциального уравнения второго порядка

Здесь v - порядок функции. Функции Бесселя могут быть первого J ±„ (х), второго y±v {х) и третьего ±v (х) родов. Они связаны следующими соотношениями:

П (х) = (/. (х) cos (vn)-/-v (x))/sm (vsi),

(6.3)

Hi (x)=y. {x)±jY.{x).

Соотношение (6.3) справедливо при дробном порядке V. При целых v оно заменяется предельным переходом.

Решениями дифференциального уравнения

dw . 1 dw

dx хЧ

к. (X) :

(6.4)

являются модифицированные функции Бесселя /v {х) и Kv (х):

п /-у(х)-/у [х) 2 sin (vo-t)

Функции Бесселя первого рода /„ {х) и In (х) при целом порядке v = n (п = 0, 1,2,...) могут вычисляться разложениями в ряды

/- м=(т)"

(/4)-/!(/-fn)!

1 = 0

Ряды для /„ {х) и /„ (х) различаются лишь знаком у члена (х/4). Это позволяет вычислять данные функции по одной программе. Программа 6.12.

Пример. /о (0,5) = 0,9384698072,

Узо (20) = 1,240153633-10"*, /„ (2) =

= 2,279585302 и /, (2) = 1,590636855.

Функции Бесселя J„ (х) и Л. (х) при дробном порядке v вычисляются по следующим разложениям в ряд:

(.=(т)

(74)

KIV (v+K+l)

Для вычисления функции Г (v-f-K-f-!) = = Г (v)-v (v + l) (v-f2) ... (v + K) должна использоваться часть программы для вычисления гамма-фуикции, например, по аппроксимации. Вычислив функции /±» [х) и-/±v (х) по формулам (6.3) и (6.4), находим функции Бесселя второго рода Y„ (х) и

К. {X).

Программа 6.13.

Пример. Для v=V=l/3 и х=1,Ъ, задав код О, получим Г (v) =2,678938244,

(х) =0,637132706!, Г (-v) =-4,062354258, / v {х) =0,2348995028 и У, [х) =0,0966101570!. Задав код 1, получим /v (л:) = !,501429!63, / v (х) = 1,622808!2 и /Cv (х) =0,2201570763.

Вычисление функций Бесселя У» {х} и Kv (х) при целых v = n невозможно. На практике приближенные значения этих функций можно найти, задав v = «±6, где б - малая величина (6» 1 -10--=-! • 10""*). Результат можно уточнить, взяв среднее при v = = п + б и v = n -6. Следует отметить, что при vO значения х9 для функций /v (х), при этом погрешность вычислений не превышает 2-10-* (с ростом v диапазон значений х расширяется и при v = 100 x<50). Погрешность вычисления У» [х) и К„ (х) выше; например, в приведенном примере Уу {х) верны 4 цифры после запятой.

Частными решениями дифференциального уравнения

"2 X dt

о) = 0.

п = 0, ±!, ±2

05 PRINT ВЫЧИСЛЕНИЕ ФУНКЦИЙ JN<X> И IN<X>

10 INPUT ВВЕДИТЕ -1 ДЛЯ JN<X> ИЛИ 1 ДЛЯ IN<X>Z

20 INPUT ВВЕДИТЕ N=N

30 INPUT ВВЕДИТЕ Х=Х

40 LETS=0! LETA=1

50 IF Н<=1 THEN 70

60 FOP 1=1 TO N! ЬЕТА=ЙЖ1: NEXT I

70 LETV=(X>2>2:LETI=e! LETB=1: LETC=1! LETR=1

80 LET 1=1+1 !LETB=I*B! LETC=«:n+I>*C

90 LETR=Z*R! LETE=<V4>*R>B>C! LETS=S+E

100 IF ABS<E>>lE-9 THEN 80

110 LETJ=<l+S>5<i<<X.2)"-N)/A

120 IF Z>0 THEN 140

130 PRINT JN«:X>=J! C-iOTO 30

140 PRINT IN<X>=J! C-iOTO 30SENII



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 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78



0.0084
Яндекс.Метрика