La función gamma

Factorial de un número

Recordaremos en primer lugar, el concepto de factorial de un número entero positivo n.

n!=n·(n-1)·(n-2)...3·2·1

MATLAB dispone de una función denominada factorial, lo calculamos también con la función prod.

>> 1*2*3*4*5
ans =   120
>> factorial(5)
ans =   120
>> prod(1:5)
ans =   120

La función gamma

Se define del siguiente modo

Γ(z)= 0 e t t z1 dt

Donde z puede ser un número real o complejo.

Esta función es más complicada en la parte negativa de z, tiene asíntotas verticales para z=0, -1, -2, -3...

syms x;
ezplot(gamma(x),[-4,5]), 
grid on

Representamos gráficamente la función gamma para z>0, llamando a la función del mismo nombre disponible en MATLAB.

>> fplot(@gamma,[0.1,4])
>> grid on
>> xlabel('x')
>> ylabel('\gamma(x)')
>> title('Función gamma')

Calculamos el mínimo de la función aproximadamente.

>> x=1:0.01:2;
>> y=gamma(x);
>> yp=diff(y);
>> idx=find(abs(yp)<0.0001);
>> x(idx)
ans =    1.4500    1.4600
>> y(idx)
ans =    0.8857    0.8856

El mínimo se encuentra entre 1 y 2 en la posición x=1.46 y vale 0.8856.

Relaciones

Relacionamos Γ(z+1) con Γ(z)

Γ(z+1)= 0 e t t z dt

Integrando por partes

u·dv =uv v·du u= t z du=z t z1 dt dv= e t dtv= e t Γ(z+1)= 0 e t t z dt = [ e t t z ] 0 +z 0 e t t z1 dt =z 0 e t t z1 dt Γ(z+1)=zΓ(z)z>0

Mediante esta relación calculamos el valor de la función gamma, para z=1,2,3,4,...

Γ(1)= 0 e t dt =1 Γ(2)=1·Γ(1)=1 Γ(3)=2·Γ(2)=2 Γ(4)=3·Γ(2)=6 ..... Γ(n+1)=n·Γ(n)=n!

>> z=1:6;
>> gamma(z)
ans =     1     1     2     6    24   120

Calculamos la función Γ(z) para valores fraccionarios de z. Por ejemplo, para z=1/2

Γ(1/2 )= 0 e t t 1/2 dt

>> syms t;
>> y=exp(-t)/sqrt(t);
>> int(y,0,inf)
ans =pi^(1/2)

Para calcular el valor de la integral hacemos el cambio t=y2.

Γ(1/2 )=2 0 exp( y 2 )dy

Se puede consultar en la tablas el valor de esta integral, o resolverla del siguiente modo

I 2 =( 0 exp ( α x 2 )dx )( 0 exp ( α y 2 )dy )= 0 0 exp( α( x 2 + y 2 ) ) dx·dy

Cambiamos de coordenadas rectangulares a polares. En coordendas polares el elemento diferencial de área dx·dy=(r·dθ)dr

x 2 + y 2 = r 2 dx·dy=(r·dθ)dr 0 π/2 0 exp( α r 2 ) r·dr·dθ= π 2 0 exp( α r 2 )·r·dr= π 4α I= 0 exp ( α x 2 )dx= 1 2 π α

Finalmente, utilizando la relación anterior Γ(z+1)=z·Γ(z)

Γ(1/2 )= π Γ(3/2 )= 1 2 Γ(1/2 )= 1 2 π .... Γ( n+ 1 2 )= 1·3·5···(2n1) 2 n π

>> n=6;
>> num=2*(1:n)-1;
>> prod(num)*sqrt(pi)/2^n
ans =  287.8853
>> gamma(n+1/2)
ans =  287.8853

Fórmula de Euler

Γ(z)Γ(1z)= π sin( πz )

>> z=1.3;
>> gamma(z)*gamma(1-z)
ans =   -3.8832
>> pi/sin(pi*z)
ans =   -3.8832

La función beta

La función beta está relacionada con la función gamma

B(m+1,n+1)= 0 1 t m ( 1t ) n dt= Γ(m+1)Γ(n+1) Γ(m+n+2)

Periodo de las oscilaciones

Consideremos una partícula de masa m oscilando a lo largo del eje X bajo la acción de una fuerza conservativa cuya energía potencial tiene la forma

E p (x)=b | x | n ,n2

b es un número real positivo y n suele ser un entero, pero puede ser un número real positivo

Representamos la energía potencial Ep(x) entre -1.2<x<1.2, tomando b=1, para cada uno de los siguientes casos: n=1,2,3,4,5. Se trata de una función simétrica, por lo que solamente es necesario calcular los valores entre 0 y 1.2

x=linspace(0,1.2,200);
hold on
for n=1:5
    y=x.^n;
    yy=[fliplr(y),y];
    xx=[-fliplr(x),x];
    plot(xx,yy,'displayName',num2str(n))
end
hold off
legend('-DynamicLegend','location','best')
grid on
ylim([0,1.2])
xlabel('x')
ylabel('E_p(x)')
title('Energía potencial')

La energía constante E de la partícula es

E= 1 2 m ( dx dt ) 2 +b | x | n

Para calcular el periodo P de la oscilación despejamos dt e integramos

dt= dx 2 m ( Eb | x | n ) P=4 0 A dx 2 m ( Eb x n ) = 4 m 2E 0 A dx 1 b x n E ,x0

El periodo es cuatro veces el tiempo que tarda en desplazarse desde el origen hasta A, la amplitud de la oscilación, la posición en la que la velocidad dx/dt de la partícula es nula. E=bAn

En el intervalo 0≤xA, x es positivo por lo que |x|=x

Para resolver la integral hacemos el cambio de variable

x= ( E b ) 1/n sin 2/n θ dx= ( E b ) 1/n 2 n sin 2/n1 θcosθ·dθ

Por otra parte,

x n = E b sin 2 θ

cuando x=A, sin2θ=1, θ=π/2

En términos de la nueva variable θ, el periodo se expresa

P= 8 n m 2E ( E b ) 1/n 0 π/2 sin 2/n1 θ·dθ

Buscamos esta integral en una tabla de integrales, pág. 397, 3.621, n° 1

0 π/2 sin μ1 θ·dθ= 2 μ2 B ( μ 2 , μ 2 )

con μ=2/n. Donde B es la función beta, teniendo en cuenta su relación con la función gamma

B(m,n)= Γ(m)Γ(n) Γ(m+n)

Obtenemos la fórmula del periodo en términos de la energía total E de la partícula

P= 8 n m 2E ( E b ) 1/n 2 2/n2 Γ 2 ( 1 n ) Γ( 2 n ) = 2 2/n n 2m E ( E b ) 1/n Γ 2 ( 1 n ) Γ( 2 n )

En términos de la amplitud A

P= 2 2/n n 2m b A 1n/2 Γ 2 ( 1 n ) Γ( 2 n )

Representamos P m b , en función de la amplitud A

hold on
for n=1:4
    f=@(x) 2^(2/n)*sqrt(2)*x.^(1-n/2)*gamma(1/n)^2/(n*gamma(2/n));
    fplot(f,[0,5],'displayName',num2str(n))
end
hold off
grid on
legend('-DynamicLegend','location','best')
ylim([0,25])
xlabel('A')
ylabel('P')
title('Periodo en función de la amplitud')

Para n=2, el periodo es independiente de la amplitud

Ejemplos

Ejemplos en el Curso de Física

Vaciado de un depósito abierto.

Referencias

Special functions and their applications. N.N. Lebedev. Prentice Hall Inc. 1965

Anyi Amezquita, María Delgado, Diego Rasero. Sistemas oscilantes intrínsecamente no lineales. Revista Brasileira de Ensino de Física, vol. 43, e20210110 (2021)

I. S. Gradshteyn, I. M. Ryzhik. Table of Integrals, Series, and Products. Seventh Edition. Elsevier (2007). Pág. 397, 3.621, n° 1