Here is mpmath's evaluation for two Struve function values:
>>> struvel(1+j, 100j)
(0.1745249349140313153158106 + 0.08029354364282519308755724j)
>>> struvel(1+j, 700j)
(-0.1721150049480079451246076 + 0.1240770953126831093464055j)
The same values in Mathematica:
In[52]:= N[StruveL[1+I, 100I], 25]
Out[52]= 0.1745249349140313153158107 + 0.0802935436428251930875572 I
In[53]:= N[StruveL[1+I, 700I], 25]
Out[53]= -0.2056171312291138282112197 + 0.0509264284065420772723951 I
I'm almost certain that the second value returned by Mathematica is wrong. The value from mpmath agrees with a high-precision direct summation of the series defining the Struve L function, and even Mathematica gives the expected value if one rewrites the L function in terms of the H function:
In[59]:= n=1+I; z=700I
Out[59]= 700 I
In[60]:= N[-I Exp[-n Pi I/2] StruveH[n, I z], 25]
Out[60]= -0.1721150049480079451246076 + 0.1240770953126831093464055 I
Maple also agrees with mpmath:
> evalf(StruveL(1+I, 700*I), 25);
-0.1721150049480079451246076 + 0.1240770953126831093464055 I
So unless Mathematica uses some nonstandard definition of Struve functions, unannounced, this very much looks like a bug in their implementation.
Wolfram Alpha reproduces the faulty value, so this still appears to be broken in Mathematica 7.
7 comments:
Great job! :)
I would be fair to tell Wofram about these bug (if not already done ;))
Xavier
Sometimes Mathematica has a bug, sometimes it is just difficult to understand how Mathematica behaves internally.
Both things are bad but to be fair with Wolfram, this is not a bug *once* you know how Mathematica works
In[8]:= StruveL[N[1+I,1000],N[700I,1000]] //N
Out[8]= -0.172115+0.124077 I
We can claim that Mathematica has no bugs at all, but just an impossible to understand semantics: If Mathematica output is wrong then change the input. If you don't like the answer, change the question. That seems to be philosophy behind.
It's apparent that Mathematica uses one algorithm at low precision and another at high precision. So forcing it to work at high precision clearly works. But the fact that the result is wrong at low precision is a bug, no doubt about it.
first let me simplify the input
In[1]:= N[StruveL[1 + I, 700 I], 1000]
also gives the right answer.
and no, the contrary is apparent. There is no algorithm switching: if I use N[ ..., 100] (which is still high precision) the erroneous result is still there. So it is not a different algorithm but it is an algorithm that depends on very high (and explicit) precision in that range of input values.
The fact that Maple (and others) gives the correct answer with *similar* input doesn't mean that this is Mathematica bug. For this case in particular it means that the rules for automatic precision are different. And at most we can say that the internal (numerical) implementation of StruveL is not optimal (outperformed by Maple).
Mathematica doesn't warranty correct numerical precision evaluation for an arbitrary function in an arbitrary range and I doubt any program does warranties that.
What to call a bug and what not in these complicated "operating" systems it a matter of taste.
The reason you got the right answer with N[..., 100] is that Mathematica had cached the high precision result. See here:
remote3:frejohl:[~]$ math
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.
In[1]:= N[StruveL[1+I, 700I], 1000];
In[2]:= N[StruveL[1+I, 700I], 100]
Out[2]= -0.172115004948007945124607606119825374911781438402257815666327260441\
> 2284212559682974833533255929038377 +
> 0.1240770953126831093464054637173332234017937943136204027954644809601120\
> 236785100654592241257674105031 I
In[3]:= Exit
Then from scratch:
remote3:frejohl:[~]$ math
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.
In[1]:= N[StruveL[1+I, 700I], 100]
Out[1]= -0.205617131229113828211219653277451317225185342987579213798338606777\
> 2107306337265415468340083692300517 +
> 0.0509264284065420772723950655021084490042078638233752974859783195081917\
> 512281936143415836657913062552 I
Getting entirely different results depending on the order of evaluations is definitely a bug in itself.
I know that Mathematica doesn't "warrant" anything (nor does mpmath), but this is not an excuse for being wrong.
Based on my own experience implementing special functions, I strongly suspect that this is a case of Mathematica implementing the wrong asymptotic expansion (the correct expansion around 0 only being used when precision is high enough) and not just of having the internal precision slightly wrong.
> The reason you got the right answer with N[..., 100] is that ...
I didn't say that I got the right answer with N[..,100], I said that I got the *wrong* answer.
Again, N[...,100] gives the wrong answer and N[...,1000] gives the right answer and I think the algorithm is the same.
> Getting entirely different results depending on the order of evaluations is definitely a bug in itself.
I agree. I agree! I was not pointing at that nor used caching in my MMA session.
> Based on my own experience
implementing special functions, strongly suspect that this is a case of Mathematica implementing the wrong asymptotic expansion (the correct expansion around 0 only being used when precision is high enough) and not just of having the internal precision slightly wrong.
So do you thing that MMA uses two different expansion depending on the precision (e.g. different expansion if it is 100 o 1000)? That seems unlikely.
Cheers,
Alfredo
I think it uses different expansions depending on both the precision and the argument. This is common operation for special functions; there's no reason at all why it should be unlikely.
This theory is further supported by benchmarking. If you try evaluating StruveL[1+I, x I] for increasing x, it gradually gets slower and slower until it hits a point where it suddenly gets much faster and remains fast for much larger x. This is indicative of a cutoff between a convergent expansion around zero and an asymptotic expansion around infinity.
I can't make a plot right now, but the performance curve should look something like the graphs here: http://fredrik-j.blogspot.com/2009/07/improved-incomplete-gamma-and.html
Post a Comment