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:
- 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:
- function mathfunc/bessi0_s(x)
- Parameters:
x [real,in]
- Return:
bessi0_s [real]
- Called from:
- Call to:
- function mathfunc/bessi1_s(x)
- Parameters:
x [real,in]
- Return:
bessi1_s [real]
- Called from:
- Call to:
- function mathfunc/bessk0_s(x)
- Parameters:
x [real,in]
- Return:
bessk0_s [real]
- Call to:
- function mathfunc/bessk1_s(x)
- Parameters:
x [real,in]
- Return:
bessk1_s [real]
- Call to:
- 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:
- function mathfunc/poly(x, coeffs)
- Parameters:
x [real]
coeffs (5) [real]
- Return:
poly [real]
- Called from: