3
$\begingroup$

I am trying to do a plot of integrals of Mathieu functions, which are built in functions in Mathematica. My problem is that it seems that Plot gives the wrong result with a factor of 2 in my case.

My code is

Clear["Global`*"];
q[λ_, L_] := (λ*L^2)/(4*π^2)
al[l_, λ_, L_] := MathieuCharacteristicA[l, q[λ, L]]
MCos[l_, λ_, L_, z_] := 
MathieuC[al[l, λ, L], q[λ, L], z]
Dp[k_, l_, λ_, L_] := (2/(π*L))*
NIntegrate[(Exp[-I*k*z])*MCos[0, λ, L, z]*
MCos[l, λ, L, z], {z, 0, π}]

Plot[{
Sqrt[q[λ, 1]] Abs[Dp[1, 1, λ, 1]]^2,
Sqrt[q[1000, 1]] Abs[Dp[1, 1, 1000, 1]]^2},
{λ, 10, 1000}, PlotRange -> {0, 0.6}]

enter image description here

Analytically, I know that the exact limit is 1/4, and not 1/2 which the Plot command gives. What am I doing wrong here?

$\endgroup$
1
$\begingroup$

This is a long comment pointing out what fails here (but with no solution proposed, though).

The issue is in NIntegrate (the same applies to Integrate, but the computations take more time). The functions q, al and MCos seem to be ok.

Dp itself is a very computationaly expensive function, which gives not very numericaly trustworthy values:

Abs[Dp[1, 100,  1, 1]]^2
Abs[Dp[1, 100., 1, 1]]^2

7.90282*10^-9

7.90282*10^-9

Abs[Dp[1, 300,  1, 1]]^2
Abs[Dp[1, 300., 1, 1]]^2

9.75498*10^-11

9.75498*10^-11

Abs[Dp[1, 1000,  1, 1]]^2
Abs[Dp[1, 1000., 1, 1]]^2

7.90139*10^-13

no result in a reasonable time

They are not trustworthy because they're are very small and likely to lead to numerical artifacts, and

f[a_] := Sqrt[q[a, 1]] Abs[Dp[1, 1, a, 1]]^2

works very fast - much faster than Dp alone.

Now, the issue is in NIntegrate because

{ f[100],  f[100.]}
{ f[300],  f[300.]}
{f[1000], f[1000.]}

{0.255864, 0.511728}

{0.249298, 0.498596}

{0.249872, 0.499745}

Also

{f[201/2], f[100.5]}

{0.255729, 0.511459}

Exact values give proper results, numerical give erroneous ones.

ListPlot[Table[f[i], {i, 10, 300, 10}], PlotRange -> {0, 0.6}]

enter image description here

but

ListPlot[Table[f[i], {i, 10., 300, 10}], PlotRange -> {0, 0.6}]

enter image description here

and Plot uses numerical values:

Plot[f[a], {a, 10, 300}, PlotRange -> {0, 0.6}]

enter image description here

A ListPlot for other exact values is also correct:

ListPlot[Table[f[i], {i, 10, 300, 1/3}], PlotRange -> {0, 0.6}]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for the answer! I did not understand that this was the difference in the Plot command, but your examples makes it very clear. It is a bit strange that even for small values of the input parameter, where the integrand is actually well-behaved, there is still a large difference between 10 and 10.0. However, I guess that is a different question. $\endgroup$ – Daniel Nov 14 '16 at 22:32
0
$\begingroup$

It seems like it is not a problem of Plot but of NIntegrate in your definition of Dp. If you use just Integrate and evaluate the expressions at the end numerically, it works fine (but takes very long).

Clear["Global`*"];
q[λ_, L_] := (λ*L^2)/(4*π^2)
al[l_, λ_, L_] := MathieuCharacteristicA[l, q[λ, L]]
MCos[l_, λ_, L_, z_] := MathieuC[al[l, λ, L], q[λ, L], z]
Dp[k_, l_, λ_, L_] := 
  (2/(π*L))*Integrate[(Exp[-I*k*z])*MCos[0, λ, L, z]*MCos[l, λ, L, z], {z, 0, π}]
f[l_] := N[Sqrt[q[l, 1]] Abs[Dp[1, 1, l, 1]]^2]
data = Table[{x, f[x]}, {x, 0, 10^3, 10^2}];
ListPlot[data, PlotRange -> All]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for the suggestion! However, I still think the problem is with Plot. If I replace Integrate with NIntegrate in your code, I get the same plot as you, but much quicker. Thus, calling the function pointwise works, and gives the correct result. $\endgroup$ – Daniel Nov 13 '16 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.