utils_mathfunc

Module

Description

Contains utility routines for mathematical functions

References:

None

Owner:

Daniel Price

Runtime parameters:

None

Dependencies:

None

Quick access

Routines:

bessi0_s(), bessi1_s(), bessk0_s(), bessk1_s(), gegenbauer_poly(), ik01a(), legendre_associated(), poly()

Variables

Subroutines and functions

subroutine  mathfunc/gegenbauer_poly(n, alpha, x, cx)

GEGENBAUER_POLY computes the Gegenbauer polynomials C(I,ALPHA,X).

Taken from the Polpak http://people.sc.fsu.edu/~jburkardt/

Discussion:

The Gegenbauer polynomial can be evaluated in Mathematica with the command

GegenbauerC[n,m,x]

ALPHA must be greater than -0.5.

If ALPHA = 1, the Gegenbauer polynomials reduce to the Chebyshev polynomials of the second kind.

Differential equation:

(1-X*X) Y’’ - (2 ALPHA + 1) X Y’ + N (N + 2 ALPHA) Y = 0

Recursion:

C(0,ALPHA,X) = 1, C(1,ALPHA,X) = 2*ALPHA*X C(N,ALPHA,X) = ( (2*N-2+2*ALPHA) * X * C(N-1,ALPHA,X)

  • ( -N+2-2*ALPHA) * C(N-2,ALPHA,X) ) / N

Norm:

Integral ( -1 <= X <= 1 )

( 1 - X^2 )^( ALPHA - 0.5 ) * C(N,ALPHA,X)^2 dX

= PI * 2^( 1 - 2 * ALPHA ) * Gamma ( N + 2 * ALPHA )

/ ( N! * ( N + ALPHA ) * ( Gamma ( ALPHA ) )^2 )

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

06 August 2004

Author:

John Burkardt

Reference:

Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65.

Parameters:

Input, integer ( kind = 4 ) N, the highest order polynomial to compute. Note that polynomials 0 through N will be computed.

Input, real ( kind = 8 ) ALPHA, a parameter which is part of the definition of the Gegenbauer polynomials. It must be greater than -0.5.

Input, real ( kind = 8 ) X, the point at which the polynomials are to be evaluated.

Output, real ( kind = 8 ) CX(0:N), the values of the first N+1 Gegenbauer polynomials at the point X.

Parameters:
  • n [integer,in]

  • alpha [real,in]

  • x [real,in]

  • cx (1 + n) [real,out]

Called from:

wang_bar()

subroutine  mathfunc/legendre_associated(n, m, x, cx)

LEGENDRE_ASSOCIATED evaluates the associated Legendre functions.

Differential equation:

(1-X*X) * Y’’ - 2 * X * Y + ( N (N+1) - (M*M/(1-X*X)) * Y = 0

First terms:

M = 0 ( = Legendre polynomials of first kind P(N,X) )

P00 = 1 P10 = 1 X P20 = ( 3 X^2 - 1)/2 P30 = ( 5 X^3 - 3 X)/2 P40 = ( 35 X^4 - 30 X^2 + 3)/8 P50 = ( 63 X^5 - 70 X^3 + 15 X)/8 P60 = (231 X^6 - 315 X^4 + 105 X^2 - 5)/16 P70 = (429 X^7 - 693 X^5 + 315 X^3 - 35 X)/16

M = 1

P01 = 0 P11 = 1 * sqrt(1-X*X) P21 = 3 * sqrt(1-X*X) * X P31 = 1.5 * sqrt(1-X*X) * (5*X*X-1) P41 = 2.5 * sqrt(1-X*X) * (7*X*X*X-3*X)

M = 2

P02 = 0 P12 = 0 P22 = 3 * (1-X*X) P32 = 15 * (1-X*X) * X P42 = 7.5 * (1-X*X) * (7*X*X-1)

M = 3

P03 = 0 P13 = 0 P23 = 0 P33 = 15 * (1-X*X)^1.5 P43 = 105 * (1-X*X)^1.5 * X

M = 4

P04 = 0 P14 = 0 P24 = 0 P34 = 0 P44 = 105 * (1-X*X)^2

Recursion:

if N < M:

P(N,M) = 0

if N = M:

P(N,M) = (2*M-1)!! * (1-X*X)^(M/2) where N!! means the product of all the odd integers less than or equal to N.

if N = M+1:

P(N,M) = X*(2*M+1)*P(M,M)

if M+1 < N:

P(N,M) = ( X*(2*N-1)*P(N-1,M) - (N+M-1)*P(N-2,M) )/(N-M)

Special values:

P(N,0,X) = P(N,X), that is, for M=0, the associated Legendre function of the first kind equals the Legendre polynomial of the first kind.

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

14 May 2004

Author:

John Burkardt

Reference:

Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34.

Parameters:

Input, integer ( kind = 4 ) N, the maximum first index of the Legendre function, which must be at least 0.

Input, integer ( kind = 4 ) M, the second index of the Legendre function, which must be at least 0, and no greater than N.

Input, real ( kind = 8 ) X, the point at which the function is to be evaluated. X must satisfy -1 <= X <= 1.

Output, real ( kind = 8 ) CX(0:N), the values of the first N+1 functions.

Parameters:
  • n [integer,in]

  • m [integer,in]

  • x [real,in]

  • cx (1 + n) [real,out]

Called from:

wang_bar()

function  mathfunc/bessi0_s(x)
Parameters:

x [real,in]

Return:

bessi0_s [real]

Called from:

bessk0_s()

Call to:

poly()

function  mathfunc/bessi1_s(x)
Parameters:

x [real,in]

Return:

bessi1_s [real]

Called from:

bessk1_s()

Call to:

poly()

function  mathfunc/bessk0_s(x)
Parameters:

x [real,in]

Return:

bessk0_s [real]

Call to:

bessi0_s(), poly()

function  mathfunc/bessk1_s(x)
Parameters:

x [real,in]

Return:

bessk1_s [real]

Call to:

bessi1_s(), poly()

subroutine  mathfunc/ik01a(x, bi0, di0, bi1, di1, bk0, dk0, bk1, dk1)

implicit double precision (a-h,o-z)

Parameters:
  • x [real,in]

  • bi0 [real,out]

  • di0 [real,out]

  • bi1 [real,out]

  • di1 [real,out]

  • bk0 [real,out] :: added by DJP, fixes compiler warning

  • dk0 [real,out]

  • bk1 [real,out]

  • dk1 [real,out]

Called from:

kfdiscsp()

function  mathfunc/poly(x, coeffs)
Parameters:
  • x [real]

  • coeffs (5) [real]

Return:

poly [real]

Called from:

bessi0_s(), bessi1_s(), bessk0_s(), bessk1_s()