Library for computation of WKB approximations of some special function

Library for computation of WKB approximations of some special functions


This library implements application of the WKB method to linear systems of ordinary differential equations (ODE) . It is numerical illustrations to the theory, developed by M. V. Fedoryuk in [1], and applied to the classic ODE for special functions.
Consider the system of n equations

(in our case of special functions n = 2). According to [1] first term of formal asymptotic expansion of solution can be written as


)dt)(ej(t)+O(l-1), j = 0, 1, ...n-1
where pj(t) - eigenvalues, determined by
and ej* - left, ej - right eigen vectors
(x))ej(x)=0, ej*(x)(
The point x0 is called a turning point of system (1.1) if equation (1.2) has a multiple root at x=x0. In the vicinity of turning point expansion (1.2) fails. At this stage we provide results with vicinities of turning points excluded, so we build (1.2) at the both sides of turning point. Points xs where det[^A](xs)=0 - singular points, are also excluded. Asymptotic expansions valid at turning points are called uniform expansions and we plan to code them later.
Recent results are implemented in static library written on Microsoft Visual C++ v6, with the header file wkblib.h displayed below:
#pragma once

#define WKB_BESSEL                    129
#define WKB_LEGENDRE                  131
#define WKB_KELVIN                    132
#define WKB_SPHERICAL_BESSEL          133
#define WKB_RICCAT_IBESSEL            135
#define WKB_AIRY                      136
#define WKB_KUMMER                    137
#define WKB_WHITTAKER                 138
#define WKB_HYPERGEOMETRIC            139
#define WKB_COULOMB                   140
#define WKB_PARABOLIC_CYLINDER        141
#define WKB_SPHEROIDAL0               143
#define WKB_SPHEROIDAL1               1431
#include <complex>
using namespace std;

typedef struct Complex2DVector
	complex<double> f;
	complex<double> g;

typedef struct
	int I;
	double K;
	complex<double> Mu;
	complex<double> Nu;
	complex<double> A;
	complex<double> B;
	complex<double> C;

void InitWKBSolutions(
				  short t,              //[in] Special function type(WKB_BESSEL, WKB_KUMMER,..)  
				  double w0,            //[in] Start point 
				  double w1,            //[in] End point
				  int N,                //[in] Number of points
				  SpecFunC SF,          //[in] Parameters of special functions  
				  int *map_size,        //[out] number of intervals ends, where solution is computed
				  double *map,          //[out] array of intervals where solution is computed
				  int *n_map            //[out] array of numbers of points on each intervals 
void WKBSolutions(
				  CV2*** wkb            //[out] 2 WKB solutions with derivatives
Library includes two functions InitWKBSolutions() and WKBSolutions().
First function tunes the code for particulat task and second performs the computations. Arguments of InitWKBSolutions() function are:
  1. short t - [input] -Special function type, one of the constants defined at the begining of header file, i.e. WKB_BESSEL, WKB_AIRY etc.
  2. double w0, w1 - [input] -left and right sides of the interval, where solution is looked for
  3. int N - [input] -number of argument values at the interwal [w0, w1]
  4. SpecFunC SF - [input] -structure , defined in wkblib.h . Different fields of this structure keep constant coefficients for different special functions, as it specified in Table 1 below.
    Function I K m n A B C
    Bessel +
    Spherical Bessel +
    Modified Spherical Bessel +
    Riccati Bessel +
    Legendre + +
    Kelvin +
    Kummer + +
    Whittaker + +
    Hypergeometric + + +
    Coulomb + +
    Parabolic cylinder +
    Spheroidal + + +
    Spheroidal + + +
    Table 1: Coefficients for special functions (see details in a description of particular special function below)
  5. int *map_size - [output] - pointer to the variable containing ,,map size''- number of ends of sub intervals of original interval [w0, w1], where solution is computed. If initial interval do not contain turning points and singular points, them *map_size = 2, or ,,map'' contain 2 points w0, w1. If there is one of those points on the interval, *map_size = 4, or initial interval id divided on two different intervals, each with two ends, with small vicinity, containig ,,bad'' point between them.
  6. *map - [output] - array of left and right ends of sub intervals, generated by turnig points and singularities
  7. *n_map - [output] - array of numbers of points of independet variable on each sub interval.
Function WKBSolutions() performs computatons. It has single argument, 3D array of 2D complex vectors wkb[m][d][n] where
0 n < map_size

, 0 d < 2, 0 n < n_map[m]
First index corresponds to different sub intervals, second index gives 2 lineary independent solutions Y0,1(x) and last index is the change of argument x. A sample program below illustrates use of the library.
#include "StdAfx.h"
#include <iostream>
#include "wkblib.h"
using namespace std;
double map_wkb[16];
int n_wkb[16];
int main(int argc, char* argv[])
	int m, n, d;
	double x0 = -0.9, x1 = 0.9, dx;
	SpecFunC sf;
	CV2 ***Solution;
	int N = 32, map_size;
	sf.Nu = complex<double>(1.0, 0.0);
	dx = (x1 - x0) / (N - 1);
	InitWKBSolutions(WKB_LEGENDRE, x0, x1, N, sf, &map_size, map_wkb, n_wkb);
	Solution = (CV2 ***)calloc(map_size / 2, sizeof(CV2**)); 
	for(m = 0; m < map_size / 2; m++)
		Solution[m] = (CV2 **)calloc(2, sizeof(CV2*));
		for(d = 0; d < 2; d++)
			Solution[m][d] = (CV2 *)calloc(n_wkb[2 * m], sizeof(CV2));
	printf("WKB Solutions\n");
	for(m = 0; m < map_size / 2; m++)
		printf("Interval %d [%E, %E], Number of Points %d\n", m, map_wkb[2 * m], map_wkb[2 * m + 1], n_wkb[2 * m]);  
		printf("     x              ReY0(x)        ImY0(x)        RedY0(x)/dx    ImdY0(x)/dx\n");
		printf("                    ReY1(x)        ImY1(x)        RedY1(x)/dx    ImdY1(x)/dx \n");
		for(n = 0; n < n_wkb[2 * m]; n++)
			printf(" %+E   %+E %+E %+E %+E\n", 
			map_wkb[2 * m] + n * dx, 
			Solution[m][0][n].f.real(), Solution[m][0][n].f.imag(),
			Solution[m][0][n].g.real(), Solution[m][0][n].g.imag()); 
			printf("%s %+E %+E %+E %+E\n","                 ", 
			Solution[m][1][n].f.real(), Solution[m][1][n].f.imag(),
			Solution[m][1][n].g.real(), Solution[m][1][n].g.imag()); 
	for(m = 0; m < map_size / 2; m++)
		for(d = 0; d < 2; d++)
	return 0;
Program output:
WKB Solutions
Interval 0 [-9.000000E-001, -8.419355E-001], Number of Points 2
     x              ReY0(x)        ImY0(x)        RedY0(x)/dx    ImdY0(x)/dx
                    ReY1(x)        ImY1(x)        RedY1(x)/dx    ImdY1(x)/dx
 -9.000000E-001   +1.000000E+000 +0.000000E+000 -1.285559E+000 +0.000000E+000
                  +1.000000E+000 +0.000000E+000 -8.188126E+000 +0.000000E+000
 -8.419355E-001   +2.631070E-002 +5.770589E-002 -4.393550E-002 -9.636143E-002
                  +3.179563E-002 -6.973570E-002 -1.307993E-001 +2.868753E-001
Interval 1 [-7.838710E-001, 7.838710E-001], Number of Points 28
     x              ReY0(x)        ImY0(x)        RedY0(x)/dx    ImdY0(x)/dx
                    ReY1(x)        ImY1(x)        RedY1(x)/dx    ImdY1(x)/dx
 -7.838710E-001   -4.697258E-002 +3.170413E-002 +6.295662E-002 -1.126780E-001
                  -4.324567E-002 -2.918866E-002 +5.796150E-002 +1.037379E-001
 -7.258065E-001   -4.855954E-002 +1.406647E-002 +5.522531E-002 -8.804912E-002
                  -4.470672E-002 -1.295040E-002 +5.084361E-002 +8.106310E-002
 -6.677419E-001   -4.701745E-002 +7.919070E-004 +5.549517E-002 -7.001045E-002
                  -4.328698E-002 -7.290753E-004 +5.109206E-002 +6.445566E-002
 -6.096774E-001   -4.359459E-002 -9.719459E-003 +5.685487E-002 -5.583833E-002
                  -4.013570E-002 +8.948295E-003 +5.234388E-002 +5.140799E-002
 -5.516129E-001   -3.897987E-002 -1.808059E-002 +5.800290E-002 -4.408337E-002
                  -3.588712E-002 +1.664603E-002 +5.340082E-002 +4.058570E-002
 -4.935484E-001   -3.364580E-002 -2.463270E-002 +5.864155E-002 -3.403970E-002
                  -3.097627E-002 +2.267828E-002 +5.398880E-002 +3.133891E-002
 -4.354839E-001   -2.794461E-002 -2.962736E-002 +5.875427E-002 -2.533097E-002
                  -2.572742E-002 +2.727666E-002 +5.409258E-002 +2.332116E-002
 -3.774194E-001   -2.214378E-002 -3.328378E-002 +5.841881E-002 -1.773375E-002
                  -2.038684E-002 +3.064298E-002 +5.378373E-002 +1.632671E-002
 -3.193548E-001   -1.644299E-002 -3.580657E-002 +5.774448E-002 -1.109829E-002
                  -1.513837E-002 +3.296560E-002 +5.316290E-002 +1.021773E-002
 -2.612903E-001   -1.098501E-002 -3.738936E-002 +5.684730E-002 -5.311346E-003
                  -1.011344E-002 +3.442280E-002 +5.233691E-002 +4.889932E-003
 -2.032258E-001   -5.864597E-003 -3.821334E-002 +5.583912E-002 -2.783396E-004
                  -5.399287E-003 +3.518141E-002 +5.140872E-002 +2.562556E-004
 -1.451613E-001   -1.136569E-003 -3.844464E-002 +5.482300E-002 +4.085007E-003
                  -1.046391E-003 +3.539436E-002 +5.047322E-002 -3.760893E-003
 -8.709677E-002   +3.176632E-003 -3.823194E-002 +5.389178E-002 +7.856280E-003
                  +2.924591E-003 +3.519854E-002 +4.961588E-002 -7.232945E-003
 -2.903226E-002   +7.078142E-003 -3.770507E-002 +5.312859E-002 +1.110770E-002
                  +6.516546E-003 +3.471347E-002 +4.891325E-002 -1.022639E-002
 +2.903226E-002   +1.059079E-002 -3.697455E-002 +5.260866E-002 +1.390643E-002
                  +9.750492E-003 +3.404091E-002 +4.843457E-002 -1.280306E-002
 +8.709677E-002   +1.375238E-002 -3.613205E-002 +5.240218E-002 +1.631461E-002
                  +1.266123E-002 +3.326525E-002 +4.824448E-002 -1.502017E-002
 +1.451613E-001   +1.661232E-002 -3.525146E-002 +5.257828E-002 +1.838945E-002
                  +1.529426E-002 +3.245453E-002 +4.840660E-002 -1.693039E-002
 +2.032258E-001   +1.922980E-002 -3.439037E-002 +5.321035E-002 +2.018379E-002
                  +1.770406E-002 +3.166176E-002 +4.898852E-002 -1.858236E-002
 +2.612903E-001   +2.167387E-002 -3.359168E-002 +5.438370E-002 +2.174716E-002
                  +1.995421E-002 +3.092644E-002 +5.006877E-002 -2.002169E-002
 +3.193548E-001   +2.402599E-002 -3.288516E-002 +5.620697E-002 +2.312779E-002
                  +2.211972E-002 +3.027597E-002 +5.174738E-002 -2.129278E-002
 +3.774194E-001   +2.638650E-002 -3.228873E-002 +5.883094E-002 +2.437610E-002
                  +2.429294E-002 +2.972687E-002 +5.416316E-002 -2.244205E-002
 +4.354839E-001   +2.888791E-002 -3.180922E-002 +6.248213E-002 +2.555098E-002
                  +2.659588E-002 +2.928541E-002 +5.752466E-002 -2.352371E-002
 +4.935484E-001   +3.172283E-002 -3.144229E-002 +6.752908E-002 +2.673244E-002
                  +2.920587E-002 +2.894759E-002 +6.217117E-002 -2.461143E-002
 +5.516129E-001   +3.520842E-002 -3.117106E-002 +7.463123E-002 +2.805218E-002
                  +3.241491E-002 +2.869788E-002 +6.870982E-002 -2.582646E-002
 +6.096774E-001   +3.996411E-002 -3.096392E-002 +8.513916E-002 +2.978778E-002
                  +3.679327E-002 +2.850717E-002 +7.838403E-002 -2.742435E-002
 +6.677419E-001   +4.758061E-002 -3.078410E-002 +1.025506E-001 +3.278698E-002
                  +4.380546E-002 +2.834162E-002 +9.441399E-002 -3.018559E-002
 +7.258065E-001   +6.575174E-002 -3.092461E-002 +1.431838E-001 +4.257600E-002
                  +6.053485E-002 +2.847099E-002 +1.318233E-001 -3.919792E-002
 +7.838710E-001   +1.000000E+000 +0.000000E+000 +2.033144E+000 +1.026534E+000
                  +1.000000E+000 +0.000000E+000 +2.033144E+000 -1.026534E+000
Interval 2 [8.419355E-001, 9.000000E-001], Number of Points 2
     x              ReY0(x)        ImY0(x)        RedY0(x)/dx    ImdY0(x)/dx
                    ReY1(x)        ImY1(x)        RedY1(x)/dx    ImdY1(x)/dx
 +8.419355E-001   +1.000000E+000 +0.000000E+000 +4.113752E+000 +0.000000E+000
                  +9.888650E-001 +0.000000E+000 +1.651278E+000 +0.000000E+000
 +9.000000E-001   +7.468092E-001 +0.000000E+000 +6.114968E+000 +0.000000E+000
                  +1.000000E+000 +0.000000E+000 +1.285559E+000 +0.000000E+000


1.Mikhail V. Fedoryuk, Asymptotic Analysis, Springer Verlag,1991

2.Handbook of Mathematical Functions, M.Abramovitz, I.A.Stegun,National Bureau of Standarda, 1964.