Expression evaluates numerically inside of Plot but not otherwise

Question

So I'm facing the somewhat bizarre situation where I can plot a function and it shows up fine, but trying to numerically integrate it causes problems. This is the plot command.

Plot[MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)],
 1/( 4 0.0275272807115071^2), x] MathieuCPrime[MathieuCharacteristicA[0,
 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x], {x, 0, 6.28},
 PlotRange -> Full]

This is what the plot looks like. But when I try to numerically integrate the same function by just replacing 'plot' above with 'NIntegrate' I get this error message:

NIntegrate::inumr: "The integrand MathieuCPrime[-623.771,329.924,x]\ MathieuS[-552.13,329.924,x] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,6.28}}."

Taking the above expression and substituting a value for x doesn't produce a number:

expr = 
  MathieuS[MathieuCharacteristicB[2 0 + 2, 
     1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), 
    x] MathieuCPrime[
    MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
    1/(4 0.0275272807115071^2), x];

expr /. x -> 0.1
(* 6.98552*10^-13 MathieuS[-552.13, 329.924, 0.1] *)

Any clues how to proceed?

Answer 1

This is not a full answer, just a long comment to give more insight into the difference between evaluating the expression within Plot or NIntegrate.

(* note the use of SetDelayed *)
expr := MathieuS[
   MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x] MathieuCPrime[
   MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x]

Substitute in a value for x first, evaluate afterwards:

Block[{x = 0.1}, expr]
(* -2.09825*10^-25 *)

Evaluate symbolically first, substitute value afterwards:

Block[{x = 0.1}, Evaluate[expr]]
(* 6.98552*10^-13 MathieuS[-552.13, 329.924, 0.1] *)

Not being familiar with this function, I do not know if there is any reason, precision related or other, which makes evaluation difficult for these particular numerical values. If there are precision related reasons, that could be important for whether the result is trustworthy ...

But this observation above does give us a clear workaround for the problem. Just prevent symbolic evaluation within NIntegrate.

f[x_?NumericQ] := 
 MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)],
    1/(4 0.0275272807115071^2), x] MathieuCPrime[
   MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x]

NIntegrate[f[x], {x, 0, 6.28}]
(* -13.2952 *)

Answer 2

Use exact numbers in expr by using Rationalize

expr = MathieuS[
    MathieuCharacteristicB[2 0 + 2, 
     1/(4 Rationalize[0.0275272807115071, 0]^2)], 
    1/(4 Rationalize[0.0275272807115071, 0]^2), x] MathieuCPrime[
    MathieuCharacteristicA[0, 1/(4 Rationalize[0.0275272807115071, 0]^2)], 
    1/(4 Rationalize[0.0275272807115071, 0]^2), x];

Plot[expr, {x, 0, 2 Pi}, PlotRange -> Full]

expr /. x -> 0.1

(*  -2.0993*10^-25  *)

NIntegrate[expr, {x, 0, 2 Pi}]

(*  -13.2952  *)

Answer 3

I don't know why this happens but a short and ugly hack would be to first evaluate the function on many points of the interval and then integrate the list of evaluated points. This can be automated with FunctionInterpolation:

f[x_] := MathieuS[
   MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x] MathieuCPrime[
   MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x]
fi = FunctionInterpolation[f[x], {x, 0, 6.28}]
NIntegrate[fi[x], {x, 0, 6.27}]

Note that FunctionInterpolation shouts out some warning messages too, but it works nonetheless.

** EDIT **

Just to add something, regarding Feyre's answer: Using Riemann sums is of course the straightforward way to go. The problem is that sampling the function at regular intervals may be oversampling it in "non-interesting" regions and undersampling in "interesting" ones. This is why I suggested to use FunctionInterpolation, as it automatically samples more points in more interesting regions.

To see this, here's a plot of the function with the interpolated grid also the density of sampled points (i.e. $1/\Delta x$):

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
xs = First@InterpolatingFunctionCoordinates[fi];
ys = InterpolatingFunctionValuesOnGrid[fi];
Plot[fi[x], {x, 0, 6.28}, 
 Epilog -> {Red, Point@Transpose@{xs, ys}}, PlotRange -> All]
ListPlot[Transpose[{MovingAverage[xs, 2], 1/Differences[xs]}]]

Answer 4

One way around this is to use Riemann Sums explicitly,

With[{Δx = 1/100}, 
 Total[Table[Δx (MathieuS[
       MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 
       1/(4 0.0275272807115071^2), x]  MathieuCPrime[
       MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
       1/(4 0.0275272807115071^2), x]), {x, 0, 
    2 Pi, Δx}]]]

-13.2952

Answer 5

Also, the problem can be solved simply by adding the pattern _Real in the function definition:

f[x_Real] := 
 MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)],
    1/(4 0.0275272807115071^2), x] MathieuCPrime[
   MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 
   1/(4 0.0275272807115071^2), x]
NIntegrate[f[x], {x, 0, 6.27}]
(*-13.2952*)