#
# Chapter 15.  Special Functions
#
#
#
# MAPLE SESSION 1
#
>    ?index,function
#
# 15.1   Overview of mathematical functions
#
#
#
# MAPLE SESSION 2
#
>    ?inifcn
#
# 15.2   Bessel functions
#
#
#
# MAPLE SESSION 3
#
>    ?Bessel
#
#
# MAPLE SESSION 4
#
>    diff(BesselJ(n,x),x);
>    evalf(BesselJ(3,0.8));
>    series(BesselJ(3,x),x,10);
#
#
# MAPLE SESSION 5
#
>    BesselJZeros(3,1);
>    evalf(%);
#
#
# MAPLE SESSION 6
#
>    zpts3 := [seq([evalf(BesselJZeros(3,n)),0],n=1..10)]: 
>    zpts4 := [seq([evalf(BesselJZeros(4,n)),0],n=1..10)]: 
>    p1 := plot(zpts3,style=point,symbol=circle,color=blue): 
>    p2 := plot(zpts4,style=point,symbol=cross,color=red): 
>    plots[display](p1,p2,axes=none);
#
# 15.3   The Gamma function
#
#
#
# MAPLE SESSION 7
#
>    plot(GAMMA(x),x=-4..4,y=-10..10,discont=true);
#
#
# MAPLE SESSION 8
#
>    GAMMA(10) = 9!;
>    GAMMA(1/2);
>    GAMMA(1/3);
>    evalf(%);
>    GAMMA(-1+I);
>    evalf(%);
#
#
# MAPLE SESSION 9
#
>    fn := n!;
>    bnk:=binomial(n,k);
>    gfn := convert(fn,GAMMA);
>    gbnk := convert(bnk,GAMMA);
#
#
# MAPLE SESSION 10
#
>    convert(gfn,factorial); 
>    convert(gbnk,binomial);
#
#
# MAPLE SESSION 11
#
>    Beta(x,y);
>    convert(%,GAMMA);
>    int(t^(x-1)*(1-t)^(y-1),t=0..1);
#
#
# MAPLE SESSION 12
#
>    lnGAMMA(-3/2);
>    evalf(%);
>    taylor(lnGAMMA(z+1),z=0,6);
#
#
# MAPLE SESSION 13
#
>    diff(lnGAMMA(z),z);
#
#
# MAPLE SESSION 14
#
>    Psi(1/6);
#
#
# MAPLE SESSION 15
#
>    Psi(1/2); 
>    Psi(1/4); 
>    Psi(3/4); 
>    exp(Psi(1/6)+Psi(5/6)+2*gamma); 
>    simplify(%);
#
#
# MAPLE SESSION 16
#
>    Psi(1,2);
#
# 15.4   Hypergeometric functions
#
#
#
# MAPLE SESSION 17
#
>    simplify( hypergeom([ ],[ ],z) );
>    simplify( hypergeom([a],[ ],z) );
>    simplify( hypergeom([ ],[3/2],-z^2/4) );
>    simplify( hypergeom([ ],[1/2],-z^2/4) );
>    simplify( hypergeom([ ],[a],-z^2/4) );
>    simplify( hypergeom([1,1],[2],-z) );
#
#
# MAPLE SESSION 18
#
>    simplify( hypergeom([a,b],[c],1) );
#
#
# MAPLE SESSION 19
#
>    simplify( hypergeom([a,b,c],[1+a-b,1+a-c],1) ); 
>    simplify( hypergeom([a,1+a/2,b,c,d],
>  [a/2,1+a-b,1+a-c,1+a-d],1) );
#
#
# MAPLE SESSION 20
#
>    Sum(pochhammer(a,n)*pochhammer(b,n)/pochhammer(c,n)/n!,
>  n=0..infinity);
>    convert(%,hypergeom);
>    simplify(%);
#
# 15.5   Elliptic integrals
#
#
#
# MAPLE SESSION 21
#
>    taylor(EllipticK(k),k=0,10);
#
#
# MAPLE SESSION 22
#
>    taylor(EllipticK(k)-Pi/2*hypergeom([1/2,1/2],[1],k^2),
>  k=0,100);
#
#
# MAPLE SESSION 23
#
>    diff(EllipticE(k),k);
#
#
# MAPLE SESSION 24
#
>    diff(EllipticK(k),k);
#
#
# MAPLE SESSION 25
#
>    k := 0.35;
>    q := EllipticNome(k);
>    EllipticModulus(q);
#
# 15.6   The AGM
#
#
#
# MAPLE SESSION 26
#
>    g1 := 1/GaussAGM(1,sqrt(2));
>    g2 := 2/Pi*int(1/sqrt(1-t^4),t=0..1);
>    evalf(g1,12);
>    evalf(g2,12);
#
#
# MAPLE SESSION 27
#
>    evalf(g1-g2,100);
#
# 15.7   Jacobi's theta functions
#
#
#
# MAPLE SESSION 28
#
>    Digits:=30: 
>    z:=rand(0..2000)()/2000.;
>    q:=rand(0..1000)()/1001.;
>    f1 := JacobiTheta1(z,q)^2*JacobiTheta4(0,q)^2: 
>    f2 := JacobiTheta3(z,q)^2*JacobiTheta2(0,q)^2: 
>    f3 := JacobiTheta2(z,q)^2*JacobiTheta3(0,q)^2: 
>    f1 - f2 + f3;
#
# 15.8   Elliptic functions
#
#
#
# MAPLE SESSION 29
#
>    u:=1.0: 
>    k:=0.5: 
>    s:=JacobiSN(u,k);
>    EllipticF(s,k);
#
#
# MAPLE SESSION 30
#
>    series(WeierstrassP(z,g2,g3),z,10);
#
#
# MAPLE SESSION 31
#
>    diff(WeierstrassP(z,g2,g3),z);
#
#
# MAPLE SESSION 32
#
>    WeierstrassPPrime(z,g2,g3)^2 - 4*WeierstrassP(z,g2,g3)^3 
>  + g2*WeierstrassP(z,g2,g3):
>    series(%,z,20): 
>    normal(%);
#
#
# MAPLE SESSION 33
#
>    diff(WeierstrassZeta(z,g2,g3),z);
#
#
# MAPLE SESSION 34
#
>    diff(WeierstrassSigma(z,g2,g3),z);
#
# 15.9   The Riemann zeta-function
#
#
#
# MAPLE SESSION 35
#
>    for n from 2 by 2 to 10 do 
>    print(zeta(n)=Zeta(n)); 
>    end do;
#
#
# MAPLE SESSION 36
#
>    plot(abs(Zeta(1/2+t*I)), t= 0..40);
#
# 15.10   Orthogonal polynomials
#
#
#
# MAPLE SESSION 37
#
>    with(orthopoly);
#
#
# MAPLE SESSION 38
#
>    with(orthopoly): 
>    P(2,a,b,x);
#
# 15.11   Integral transforms
#
#
#
# MAPLE SESSION 39
#
>    with(inttrans);
#
# 15.11.1   Fourier transforms
#
#
#
# MAPLE SESSION 40
#
>    with(inttrans): 
>    f := exp(-t^2/2);
>    F := fourier(f,t,x);
#
#
# MAPLE SESSION 41
#
>    with(inttrans): 
>    F := sqrt(2*Pi)*exp(-x^2/2);
>    f := invfourier(F,x,t);
#
#
# MAPLE SESSION 42
#
>    with(inttrans): 
>    assume(a>0): 
>    f := exp(-a*t);
>    F := fouriercos(f,t,x);
>    Phi := fouriersin(f,t,x);
#
#
# MAPLE SESSION 43
#
>    fouriercos(F,x,t); 
>    radsimp(%); 
>    fouriersin(Phi,x,t); 
>    radsimp(%);
#
# 15.11.2   Hilbert transform
#
#
#
# MAPLE SESSION 44
#
>    with(inttrans): 
>    f := 1/(1+t^2);
>    hilbert(f,t,x);
#
#
# MAPLE SESSION 45
#
>    with(inttrans): 
>    F := -x/(1 + x^2);
>    invhilbert(F,x,t);
>    hilbert(F,x,t);
#
# 15.11.3   Mellin transform
#
#
#
# MAPLE SESSION 46
#
>    with(inttrans): 
>    mellin(sin(t),t,s);
#
#
# MAPLE SESSION 47
#
>    with(inttrans): 
>    F := GAMMA(s)*sin(Pi*s/2);
>    invmellin(F,s,t);
#
#
# MAPLE SESSION 48
#
>    with(inttrans): 
>    F := GAMMA(s);
>    invmellin(F,s,t);
#
#
# MAPLE SESSION 49
#
>    with(inttrans): 
>    F := GAMMA(s);
>    invmellin(F,s,t,0..infinity);
#
#
# MAPLE SESSION 50
#
>    with(inttrans): 
>    f := exp(-t); 
>    mellin(f,t,s);
#
# 15.12   Fast Fourier transform
#
#
#
# MAPLE SESSION 51
#
>    x := array([seq(k,k=1..8)]);
>    y := array([seq(0,k=1..8)]);
>    FFT(3,x,y);
>    print(x);
>    print(y);
#
#
# MAPLE SESSION 52
#
>    omega := exp(2*Pi*I/8); 
>    sum((k+1)*z^ k,k=0..7): 
>    p := unapply(%,z); 
>    for k from 0 to 7 do 
>   k, x[k+1]+I*y[k+1], evalf(p(omega^ k));
>   od:
#
#
# MAPLE SESSION 53
#
>    x := array([seq(k,k=1..8)]): 
>    y := array([seq(0,k=1..8)]): 
>    FFT(3,x,y): 
>    iFFT(3,x,y);
>    print(x);
>    print(y);
#
# 15.13   Asymptotic expansion
#
#
#
# MAPLE SESSION 54
#
>    z:='z':  
>    asympt(Psi(z),z,3);
#
#
# MAPLE SESSION 55
#
>    asympt(GAMMA(z),z,3);