Funciones especiales y Fortran 77: I. gamma Γ(z), factoriales n!, coeficiente binomial (n k) y beta ß(z,w)

Funciones especiales y Fortran 77: I. gamma Γ(z), factoriales n!, coeficiente binomial (nk) y beta β(z, w) Alejandro R. Urzúa Pineda Grupo de Óptica Matemática Departamento de Física Teórica Instituto de Ciencias Físicas, UNAM

Resumen El siguiente trabajo presenta un resumen de las funciones especiales: gamma, beta, y el coeficiente binomial, así como su implementación en Fortran 77. Se construye primero la función gamma, enseguida se definen los factoriales n! y ln(n!), base necesaria para definir al coeficiente binomial, por último se construye la función beta.

1.

Introducción

Las funciones especiales en realidad no tienen nada de especial, sólo son un denominativo bastante altisonante que algún autor les dió al ser objetos matemáticos no triviales que no se ven todos los días. También llamadas funciones trascendentales superiores (¿superiores a qué?) o funciones de la física matemática (aunque ocurren mayoritariamente en esta área, es posible que ocurran en muchos otros lados de la ciencia). Aunque simplemente debería denominárseles como “funciones útiles"(muy útiles); son funciones mediante las cuales se puede resolver una ecuación diferencial, la búsqueda de una raíz polinomial o más aún, el tratamiento de un problema en el plano complejo, son funciones las cuales llevan un nombre único y se comportan de manera única. Dos de estas, gamma y coeficiente binomial, sirven para definir de manera recursiva una tercera función especial llamada polinomio ortogonal de Kravchuk, otro tipo de funciones que no destacaré ahora. ´ En general, este tipo de funciones hoy día se pueden evaluar numerica y simbólicamente de forma fácil en programas sofisticados de cómputo científico; pero, al ser exigido un cálculo puramente numérico con miras a la implementación en paralelo en computadoras no del tipo super-, se decidió recurrir a la implementación en el lenguaje Fortran 77, lenguaje diseñado puramente para las tareas de cálculo numérico intrínseco.

2.

Teoría de funciones e implementación computacional

2.1.

Función gamma Γ(z)

La función gamma se define por la integral ∞

Γ(z) = ∫

tz−1 e−t dt,

0

z∈Z

La función gamma se relaciona con el factorial desplazado por un paso, es decir n! = Γ(n + 1).

1

(1)

Además, cumple con la relación de recurrencia Γ(z + 1) = zΓ(z).

(2)

Si la función es conocida para z > 1 o, más generalmente, para el semiplano complejo Re(z) > 1, es posible obtener la evaluación para z < 1 o Re(z) < 1 por la fórmula de reflexión π πz Γ(1 − z) = = . (3) Γ(z) sin(πz) Γ(1 + z) sin(πz) Es de notarse que Γ(z) tiene un polo en z = 0, y en todos los valores enteros negativos de z. Existe una amplia variedad de métodos para calcular la función Γ(z) numericamente, pero ninguno es tan concreto como la aproximación derivada por Lanczos [1]. Dejando de lado los detalles técnicos de la derivación de la aproximación y, si nos centramos en el resultado: Para cierta elección de enteros γ y N , y para ciertos coeficientes cn , n = 1, ⋯, N , la función gamma está dada por la serie 1

1 z+ 2 1 Γ(z + 1) = (z + γ + ) exp (− {z + γ + }) 2 2 √ c1 c2 cN + +⋯+ + ) , (z > 0) (4) × 2π (c0 + z+1 z+2 z+N Se puede apreciar que es un tipo de aproximación de Stirling para funciones analíticas, pero con una serie de correcciones que toman en cuenta los primeros polos en el semiplano complejo derecho de gamma. La constante c0 es muy cercana a 1. El término de error es parametrizado por . Por ejemplo, para γ = 5, N = 6, con un conjunto apropiado de valores de c’s, el error es menor que ∣∣ < 2 × 10−10 . ¡Impresionante! Ahora, es más conveniente en implementar LnΓ(x) que Γ(x), esto debido a que la última forma puede ocasionar un bucle infinito en la representación de punto flotante de algunos núcleos computacionales para valores pequeños de x. Con todo esto en cuenta, tomando la Ecuación (4), es posible computar el logaritmo de gamma con sólo una llamada a la rutina propia y unas veinticinto operaciones aritméticas de bajo-medio nivel. 2.1.1.

Código de la función gamma function gammln ( x x ) integer j real gammln , x x double precision s e r , s t p , tmp , x , y , c o f ( 6 ) save cof , st p data c o f , s t p / 7 6 . 1 8 0 0 9 1 7 2 9 4 7 1 4 6 d0 , − 8 6 . 5 0 5 3 2 0 3 2 9 4 1 6 7 7 d0 , * 2 4 . 0 1 4 0 9 8 2 4 0 8 3 0 9 1 d0 , − 1 . 2 3 1 7 3 9 5 7 2 4 5 0 1 5 5 d0 , . 1 2 0 8 6 5 0 9 7 3 8 6 6 1 7 9 d − 2 , * − . 5 3 9 5 2 3 9 3 8 4 9 5 3 d − 5 , 2 . 5 0 6 6 2 8 2 7 4 6 3 1 0 0 0 5 d0 /

11

x=xx y=x tmp= x + 5 . 5 d0 tmp = ( x + 0 . 5 d0 ) * l o g ( tmp ) − tmp s e r = 1 . 0 0 0 0 0 0 0 0 0 0 1 9 0 0 1 5 d0 do 1 1 j = 1 , 6 y = y + 1 . d0 ser=ser + cof ( j ) / y enddo gammln=tmp + l o g ( s t p * s e r / x ) return end

2

2.2.

Factorial n!

Considerando al factorial n! como un subproducto del computo de la función gamma, es posible obtenerlo para valor enteros pequeños, porque, de todas forma, para valores exageradamente grandes, se tendría un bucle infinito de cálculo. 2.2.1.

Código del factorial function f a c t r l ( n ) integer n real f a c t r l integer j , n t o p real a ( 3 3 ) , gammln s a v e ntop , a data ntop , a ( 1 ) / 0 , 1 . /

10

i f ( n . l t . 0 ) then pause ’ f a c t o r i a l n e g a t i v o en f a c t r l ’ else i f ( n . l e . n t o p ) then f a c t r l =a ( n + 1 ) else i f ( n . l e . 3 2 ) then do 1 0 j = n t o p + 1 , n a ( j +1)= j *a ( j ) enddo n t o p =n f a c t r l =a ( n + 1 ) else f a c t r l = exp ( gammln ( n + 1 . ) ) endif end

2.3.

Coeficiente binomial (nk)

Una de las ventajas del cómputo del factorial n! es que es exacto para valores pequeôsos de n, ya que las multiplicaciones de punto flotante son exactas en todas las computadoras hoy día. Para los coeficientes binomiales es posible obtener cierto grado de exactitud debido a que los factoriales individuales pueden ser aislados antes de entrar en un bucle infinito. El coeficiente binomial se define por n n! ( )= , k k!(n − k)! 2.3.1.

0≤k≤n

Código del coeficiente binomial function binom ( n , k ) integer n , k real binom , f a c t l n binom = n i n t ( exp ( f a c t l n ( n ) − f a c t l n ( k ) − f a c t l n ( n−k ) ) ) return

3

(5)

end function f a c t l n ( n ) integer n real f a c t l n , a ( 1 0 0 ) , gammln save a data a / 1 0 0 * − 1 . / i f ( n . l t . 0 ) pause ’ F a c t o r i a l n e g a t i v o e n c o n t r a d o ’ i f ( n . l e . 9 9 ) then i f ( a ( n + 1 ) . l t . 0 . ) a ( n + 1 ) = gammln ( n + 1 . ) f a c t l n =a ( n + 1 ) else f a c t l n =gammln ( n + 1 . ) endif return end

Si el problema en específico requiere de series relacionadas con el coeficiente binomial, una buena idea es usar relaciones de recurrencia, por ejemplo n+1 n+1 n n n ( )= ( )=( )+( ) k n−k+1 k k k−1 n n−k n ( ) ( )= k+1 k k+1

2.4.

(6)

Función beta

Finalmente, habiendo construido las implementaciones para el cálculo de funciones combinatorias de entradas enteras, ahora es posible establecer el cálculo de la función beta, la cual se define por 1

B(z, w) = B(w, z) = ∫

tz−1 (1 − t)w−1 dt.

0

(7)

la cual está relacionada a la función gamma de la forma B(z, w) =

Γ(z)Γ(w) . Γ(z + w)

(8)

De la Ec. (8) se ve justificado el hecho de haber construido primero las funciones combinatorias y culminar con la función beta. La implementación computacional se presenta a continuación function b e t a ( z , w ) real b e t a , w , z , gammln b e t a = exp ( gammln ( z ) + gammln ( w) − gammln ( z +w ) return end

Referencias [1] Lanczos, C. SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp.86-96. (1964)

4