Prev Next Index-> contents reference index search external Up-> f2cad library dbesks f2cad-> license Install get_started prototype library run.cpp utility add2lib.sh whats_new library-> d9knus d9lgmc daxpy dbesks dbskes dcsevl ddot dgamlm dgamma dgefa dgesl dint dscal idamax initds xerror dbesks-> dbesks.f Headings-> Prototype Fortran Source Notation Description Special Case

dbesks

Prototype
int f2cad::dbesks_(doublereal *xnu, doublereal *x, integer *nin, doublereal *bk);

Fortran Source
dbesks.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 dbesks computes the values   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/dbesks.hpp> test_result dbesks(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 bk[2]; integer nin = 2; f2cad::dbesks_(&xnu, x+0, &nin, bk); // dependent variables integer m = 1; // number of dependent variables doublereal f[1]; // vector of dependent variables f[0] = bk[1]; // value f[0] = K_(nu+1) (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 ) * exp(-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 = - check; check += - term * exp(-x0) / (x0 * x0); check += - (1. + 1. / x0) * exp(-x0) * term * term * term / pi; ok &= f2cad::near_equal(p, check, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dbesks.cpp