Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

bad Mathieu functions

29 views
Skip to first unread message

becko BECKO

unread,
Mar 11, 2010, 6:35:08 AM3/11/10
to
I am no expert in Mathieu functions, but I don't think this gives the right result:

ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]

Plot[I ce1[3, q, I], {q, 0, 1000}]


In another system you get a smooth graph, making very small oscillations about zer o as q increases. I've read that Mathieu functions are difficult to deal with. I guess that Mathematica's implementation doesn't support arbitrary values of the parameters. I made a bug report in the WRI site. Maybe someone here has some comments to share? It would be nice if there were a package or something with more robust implementations of Mathieu functions.

becko BECKO

unread,
Mar 12, 2010, 7:11:49 AM3/12/10
to
Thanks Carl, Peter and Daniel.

Guess I was too hasty to send that bug report... Though it would be useful
to issue a warning or something (perhaps putting it in the Possible Issues
section of the help on Mathieu functions). I'm kind of falling in love with
Mathematica, so it's reassuring to know that this was my mistake and not
Mathematica's.

Another question. How do I know that WorkingPrecision->50 is enough? What if I wanted to Plot all the way to q=100000?


> Date: Thu, 11 Mar 2010 08:26:43 -0500
> From: ca...@wolfram.com
> To: beck...@hotmail.com
> CC: math...@smc.vnet.net
> Subject: Re: bad Mathieu functions

> The problem is that for large q ce1 will return incorrect values for
> machine numbers. On the other hand, for extended precision numbers it
> will return correct values. Compare:
>
> Machine number input:
>
> In[16]:= ce1[3, 1000., I]
>
> Out[16]= 0. - 2200.76 I
>
> Extended precision number input:
>
> In[17]:= ce1[3, 1000`50, I]
>
> Out[17]= 7.29389*10^-22 I
>
> If you want to plot this, try instead:
>
> Plot[I ce1[3, q, I], {q, 0, 1000}, PlotRange->{-5,5}, WorkingPrecision->50]
>
> Carl Woll
> Wolfram Research
>

gekko

unread,
Mar 12, 2010, 7:12:10 AM3/12/10
to


Au contraire, Mathematica is one of the few systems that does support
arbitrary precision evaluation of all special for arbitrary values of
parameters. However, there is no getting around the fact that some
regions of parameter space are just hard to compute!

In the case of Mathieu functions, the series expansion by which they
are computed converges more slowly as the parameter q becomes larger.
As a consequence, to get meaningful values for large q you need to
work with higher precision.

Compare the following three cases (using your ce1 function):

In[179]:= {ce1[3, 1000., I], ce1[3, 1000`20, I], ce1[3, 1000`60, I]}

Out[179]= {0. + 1608.56 I, 0.*10^2 I, 7.293885910335973*10^-22 I}

In order, what's going on here is:

1. Machine-precision evaluation: all internal computations are done to
machine precision, which is insufficient in this case, and hence the
final answer is (machine-precision) garbage.

2. High-precision evaluation, using 20 significant digits: here
Mathematica pays attention to significant figures and correctly
returns a zero-precision result, since 20 digits is not enough to
compute the underlying sum.

3. High-precision evaluation, using 60 significant digits: here 60
digits is high enough to compute the sum and hence return a meaningful
answer. Notice that the final answer has far fewer significant digits
than the input, since precision is lost evaluating the sum.


Now, you could legitimately argue that Mathematica should be "smarter"
in cases 1 and 2 above. For case 1, I would say that the current
behaviour reasonable (and somewhat "expected") since it is a general
design decision that machine precision computation should be fast at
the expense of complete control over precision.

For case 2 I would tend to agree that Mathematica *could* be smarter:
there is no reason why the internal algorithm couldn't recognize that
the output has insufficient precision, and thus increase the internal
working precision.


In any case, to return to your original example: if you want to get
your plot "working", you can simply force it to use higher precision:

LogPlot[Abs[I ce1[3, SetPrecision[q,100], I]], {q, 0, 1000}, PlotRange-
>All]

Hope this helps.

Cheers, P.

Carl K. Woll

unread,
Mar 12, 2010, 7:13:46 AM3/12/10
to
On 3/11/2010 6:35 AM, becko BECKO wrote:
> I am no expert in Mathieu functions, but I don't think this gives the right result:
>
> ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]
>
> Plot[I ce1[3, q, I], {q, 0, 1000}]
>
>
> In another system you get a smooth graph, making very small oscillations about zer o as q increases. I've read that Mathieu functions are difficult to deal with. I guess that Mathematica's implementation doesn't support arbitrary values of the parameters. I made a bug report in the WRI site. Maybe someone here has some comments to share? It would be nice if there were a package or something with more robust implementations of Mathieu functions.
>
>
>

The problem is that for large q ce1 will return incorrect values for

Peter Pein

unread,
Mar 12, 2010, 7:13:57 AM3/12/10
to
Am 11.03.2010 12:35, schrieb becko BECKO:
> I am no expert in Mathieu functions, but I don't think this gives the right result:
>
> ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]
>
> Plot[I ce1[3, q, I], {q, 0, 1000}]
>
>
> In another system you get a smooth graph, making very small oscillations about zer o as q increases.
...

Oh yes, these oscillations are really tiny as you can see with
Mathematica using:

Plot[I ce1[3,q,I] /. q->SetPrecision[qq,123], {qq, 750, 1000},
Frame->True, Axes->None]


DrMajorBob

unread,
Mar 13, 2010, 7:57:28 AM3/13/10
to
A fourth option is to set $MinPrecision. Help for it has the following
examples:

Block[{$MinPrecision = $MachinePrecision},
N[Pi, 50] - 265099323460521503743/84383735478118508040]

4.066722320572635*10^-41

Block[{$MinPrecision = 20, $MaxPrecision = 20},
1 - Sqrt[1 - Exp[-10`20]]]

0.000022700222531293910900

In both cases, you start with an "arbitrary precision" number, and
$MinPrecision prevents calculations from dropping to zero precision.

If you start with a machine precision number, you lose that advantage, so
that $MinPrecision has no effect in this code:

Block[{$MinPrecision = $MachinePrecision},
N@Pi - 265099323460521503743/84383735478118508040]

4.44089*10^-16

Block[{$MinPrecision = 50},
N@Pi - 265099323460521503743/84383735478118508040]

4.44089*10^-16

This throws an error, too, because you're starting with too little
precision to satisfy $MinPrecision:

Block[{$MinPrecision = 50},
N[Pi, 20] - 265099323460521503743/84383735478118508040]

4.0667223205726345135327471876400850811855821022270*10^-41

(Yet it seems to have worked.)

Bobby

On Fri, 12 Mar 2010 06:12:26 -0600, gekko <pfal...@gmail.com> wrote:

> On Mar 11, 10:35 pm, becko BECKO <becko...@hotmail.com> wrote:


--
DrMaj...@yahoo.com

DrMajorBob

unread,
Mar 13, 2010, 7:58:03 AM3/13/10
to
Or, a trifle simpler:

ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]

Plot[I ce1[3, SetPrecision[qq, 123], I], {qq, 750, 1000},


Frame -> True, Axes -> None]

or simpler yet:

Plot[I ce1[3, qq, I], {qq, 750, 1000}, Frame -> True, Axes -> None,
WorkingPrecision -> 123]

Note that this doesn't work at all, and it's much slower:

ce1[r_, q_, z_] :=

MathieuCPrime[MathieuCharacteristicA[r, SetPrecision[q, 123]], q, z]
Plot[I ce1[3, qq, I], {qq, 750, 1000}, Frame -> True, Axes -> None]

Does anyone have a clue why that is?

Bobby

On Fri, 12 Mar 2010 06:14:16 -0600, Peter Pein <pet...@dordos.net> wrote:

> Am 11.03.2010 12:35, schrieb becko BECKO:

>> I am no expert in Mathieu functions, but I don't think this gives the
>> right result:
>>
>> ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]
>>
>> Plot[I ce1[3, q, I], {q, 0, 1000}]
>>
>>
>> In another system you get a smooth graph, making very small
>> oscillations about zer o as q increases.

> ...
>
> Oh yes, these oscillations are really tiny as you can see with
> Mathematica using:
>
> Plot[I ce1[3,q,I] /. q->SetPrecision[qq,123], {qq, 750, 1000},
> Frame->True, Axes->None]
>
>


--
DrMaj...@yahoo.com

becko BECKO

unread,
Mar 13, 2010, 8:00:34 AM3/13/10
to

Nevermind. I think I'm getting the point (thanks gekko for detailed reply).
With MachinePrecision, Mathematica doesn't keep track of precision (in the interest of speed), so there is no way to issue a warning since Mathematica doesn't know the precision of the result. If I specify the precision xplicitly (using tick marks) Mathematica will keep track of precision (making it a bit slower), and then I'll know the precision of the result.


> Date: Fri, 12 Mar 2010 07:12:04 -0500
> From: beck...@hotmail.com


> Subject: Re: bad Mathieu functions

> To: math...@smc.vnet.net

David Park

unread,
Mar 13, 2010, 8:09:09 AM3/13/10
to
Use greater WorkingPrecision.

Plot[I ce1[3, q, I], {q, 0, 100},
WorkingPrecision -> 30,
PlotRange -> {-1, 1} 5]


David Park
djm...@comcast.net
http://home.comcast.net/~djmpark/

From: becko BECKO [mailto:beck...@hotmail.com]


I am no expert in Mathieu functions, but I don't think this gives the right
result:

ce1[r_, q_, z_] := MathieuCPrime[MathieuCharacteristicA[r, q], q, z]

Plot[I ce1[3, q, I], {q, 0, 1000}]


In another system you get a smooth graph, making very small oscillations

0 new messages