erf(1.0+5*%i) returns
- 2.7837671212125964e+9 %i - 1.0786562260474622e+9
But the true answer (as given by Wolfram/Alpha) is
- 2.7837702922089653e+9 %i - 1.0786931161985406e+9
Only about 4 digits are correct.
Barton Willis
2011-03-18
By the way: using a 1F1 representation for erf, the values agree with Wolfram | Alpha. But the hypergeometric code has to switch to bigfloats to do this:
(%i11) my_erf(x) := nfloat(2*x*hypergeometric([1/2],[3/2],-x^2)/sqrt(%pi),[]);
(%o11) my_erf(x):=nfloat((2*x*hypergeometric([1/2],[3/2],-x^2))/sqrt(%pi),[])
(%i12) load(hypergeometric)$
(%i14) my_erf(1.0+5*%i);
(%o14) -2.783777029220896b9*%i-1.078693116198541b9
Raymond Toy
2011-03-24
I think the inaccuracy comes from gamma-incomplete, which is used to evaluate erf. For this particular argument, gamma-incomplete uses the series.
I think the test of when to use the series could be relaxed a bit. When I set *gamma-imag* to 0.1, the continued fraction is used, which gives
- 2.7837770292209086e+9 %i - 1.0786931161985383e+9
This is better.
I've done some experiments using the continued fraction http://functions.wolfram.com/06.06.10.0005.01. This seems to work a bit better.
Dieter Kaiser
2011-03-24
Hello Ray,
you are right. The problem is the function gamma_incomplete. The continued fraction expansion of the function tends to converge to a wrong value for some ranges of the argument.
In the past, I had searched for some hints in the literature to solve this problem more general. But I had not found a simple solution. Therefore, I had done a heuristic approach and checked several ranges for the arguments to achieve a better approximation of the function gamma_incomplete.
I have to work again on this problem.
Dieter Kaiser
Raymond Toy
2011-03-24
Here is an alternative method of computing gamma_incomplete:
http://functions.wolfram.com/06.06.10.0007.01
There appear to be no restrictions on the range of validity of this continued fraction, so it can probably be used for other areas. I've tested this a bit (with oct), and it converges quickly for things on or near the negative real axis. The few tests I've done indicate that accuracy is good.
Raymond Toy
2011-03-29
Actually, the continued fraction
http://functions.wolfram.com/06.06.10.0009.01
appears to converge faster and to be more accurate. The previously mentioned fraction has some issues for points near the negative real axis where the accuracy is not as good as we might expect.
Raymond Toy
2012-08-17
Raymond Toy
2012-08-17
Appears to have been fixed some time ago.