|
Главная -> Справочник по алгоритмам 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 |
|