LabVIEW
[LabVIEW Bug Report] Parabolic Cylinder Function
I have reported on the poor implementation of the Kummer function in the past (the situation being far from resolved, despite claims to the contrary).
I am sorry to have to report on similar problems with another special function.
Introducing the parabolic cylinder function D_nu(x):
It is defined, for instance, here: http://mathworld.wolfram.com/ParabolicCylinderFunction.html
It can be conveniently calculated online using the Wolfram Alpha engine. Just type ParabolicCylinderD[nu,x] where nu is the function's parameter and x its argument.
So let's start with the usual precision game (LabvIEW 2015f2 64 bit first, then Wolfram Alpha):
- D[0,1.5] = 0.5697828240567580600000
= 0.569782824730923009766629689829122815884638474327995977291008...
difference ~O(1E-9)
- D[1.3333,1.5] = 0.8895635352834155800000
= 0.889563535907529859621528428740936371747863209252972283125817...
difference ~O(1E-9)
This is prefectly acceptable precision (considering that LV uses double precision, we might expect better, but I am not complaining at this stage).
- D[1.3333,10] = -0.111362756256787130000000
= 2.985207642444327055287247415431910447206487785233086081... × 10^-10
right... I think so too.
- D[1.3333,100] = NaN
= 8.51894064312482769698267175786164468232938918634102127... × 10^-1084
The dreaded NaN again!
In general, the limit of D_nu[x] for large x should be zero and therefore NaN is not an answer (but it is the systematic answer of NI's implementation). Probably, as for the Kummer function, the algorithm used to estimate the function is suboptimal.
As such, I would qualify the function as useless and would recommend (again) NI to put some serious effort in vetting their functions against trustworthy references.
As I am looking into a way to implement the parabolic cylinder function, I am reading Abramowitz & Stegun and, irony of ironies, it turns out that the parabolic cylinder function can be expressed in terms of... the Kummer function (which is also apparent from the Mathworld page I linked to above)! To (not) simplify things, Abramowitz & Stegun, Gradshteyn & Ryzhik (G&R), Mathworld, etc, all give different relations. I am using G&R:
D_nu[x] = 2^(nu/2)*exp(-xˆ2/4){piˆ1/2/Gamma((1-nu)/2)*M(xˆ2/2, -nu/2, 1/2) - x*(2*pi)ˆ1/2/Gamma(-nu/2)*M(xˆ2/2, (1-nu)/2, 3/2)
where M(x a, b) is the Kummer function.
Using this relation, we obtain, for the examples above (first two results as before, LV then Wolfram Alpha, 3rd line equal the right hand side of the equation above computed by LV (*), 4th computed by Mathematica with 50 digit precision (**))
- D[0,1.5] = 0.5697828240567580640000
= 0.56978282473092300976662968982912281588463847432799597729100
= 0.56978282473092289900
= 0.5697828247309230097666296898291228158846384743280
- D[1.3333, 1.5] = 0.8895635352834155850000
= 0.88956353590752985962152842874093637174786320925297228312581
= 0.88956353590752990100
= 0.88956353590752985870405653424728494877848330080
- D[1.3333,10] = -0.111362756256787127000000
= 2.985207642444327055287247415431910447206487785
= 6.7723720765602794000E-8
= 2.985207642444326573517538339703 x 10^-10
- D[1.3333,100] = NaN
= 8.5189406431248276969826717578616446823293891863
= NaN
= 0 with 50, 100 or 1000 digits of precision but with 10,000 digits:
= 8.5189406431248249373656419520216651575210966302... x 10^-1084
In other words, it is possible to get slightly better results using the imprecise Kummer function implementation of NI, but it still falls short for most practical purposes. The main reason is that the formula requires computing the Kummer function for large values (x^2/2 = 5,000), and that returns an Inf result.
(*) The VI I used for the modified Parabolic Cylinder function can be found below. Note that I had to go around the limitation of not having a Gamma function for negative values.
(**) Mathematica 9.0
As a preliminary comment, I have played a bit with some of the Trigonometric, Hyperbolic and Exponential functions, and have found them so far to return results with the expected precision (16 significant digits for DBL).
This appears not to be the case for special functions.
Take the Bessel function of the first kind for instance:
J_5(150150) = -0.000,132,928,258,3827989 (LV 2015)
= -0.000,132,928,258,390,699.65590104237064983884024569853 (Mathematica 9.0, 50 digit precision)
>>> 10 significant digits
J_50(150150) = -0.00,205,583,043,369,971,7 (LV 2015)
= -0.00,205,583,043,369,961,2701772948883589768004682846459 (Mathematica 9.0, 50 digit precision)
>> 13 significant digits
Pretty variable and consistantly degrading for large numbers...
