Prev Next

dbskes

Prototype
int f2cad::dbskes_(doublereal *xnu, doublereal *x, integer *nin, doublereal *bke);

Fortran Source
dbskes.f

Notation
We use  K_{\nu} (x) to denote the modified Bessel function of the second kind; i.e., it is the solution of the ODE  \[
x^2 y^{(2)} (x) + x y^{(1)} (x) - ( x^2 + \nu^2 ) = 0
\] 
that goes to zero as  x \rightarrow \infty .

Description
The routine dbskes computes the values  \exp(x) K_{\nu + i} (x) where  x > 0 ,  -1 < \nu < 1 ,  i = 0 , 1 , \ldots , nin-1 (if  nin > 0 ),  i = 0 , -1 , \ldots , nin+1 (if  nin < 0 ).

Special Case
By Formula 10.2.17 of Abramowitz and Stegun  \[
     K_{3/2} ( x ) = \sqrt{  \pi / ( 2 x ) } ( 1 + 1 / x ) \exp( - x )
\] 
 

# include <cmath>
# include <f2cad/dbskes.hpp>

test_result dbskes(void)
{	bool ok = true;

	// independent variable
	integer    n = 1;          // number of independent variables
	doublereal x[1];           // vector of independent variables
	x[0]         = .75;        // value of independent variable
	f2cad::Independent(n, x);  // declare x as independent variable vector

	// call the routine
	doublereal xnu = .5;
	doublereal bke[2];
	integer    nin = 2;  
	f2cad::dbskes_(&xnu, x+0, &nin, bke);

	// dependent variables 
	integer     m = 1;          // number of dependent variables
	doublereal  f[1];           // vector of dependent variables
	f[0] = bke[1];              // value f[0] = K_(nu+1) (x[0]) * exp( x[0] )
	f2cad::Dependent(m, f);     // declare f as dependent variable vector

	// check function value
	double x0    = f2cad::Value( x[0] );
	double f0    = f2cad::Value( f[0] );
	double pi    = 4. * std::atan(1.);
	double term  = 1. / sqrt( 2. * x0 / pi );
	double check = term * (1. + 1. / x0 );
	ok          &= f2cad::near_equal(f0, check, 1e-10, 1e-10);

	// Evaluate the derivative of f[0] w.r.t x[0]
	double p = f2cad::Partial<doublereal>(0, 0);
	check    = - (1. + 1. / x0) * term * term * term / pi;
	check   += - term / (x0 * x0);
	ok      &= f2cad::near_equal(p, check, 1e-10, 1e-10);

	if( ok ) 
		return test_pass;
	return test_fail;
}

Input File: example/dbskes.cpp