# Expression evaluates numerically inside of Plot but not otherwise

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?

• Nov 8 '16 at 16:32

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 *)


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  *)

• any thoughts on why this is necessary? I assume you've seen this sort of thing before since it seems like one of those tricks everyone tends to accrue. Nov 8 '16 at 16:03
• @N.J.Evans - Mathematica's internal algorithms generally operate in the complex domain. Consequently, complicated evaluations can accumulate round-off errors that result in complex artifacts. Exact numbers facilitate cancellation of the complex components in what should be real-valued expressions. Nov 8 '16 at 17:45

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["DifferentialEquationsInterpolatingFunctionAnatomy"];
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]}]]


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

• Note that the difference between this approach and mine is that FunctionInterpolation uses a "smarter" sampling of the functions, supposedly evaluating it on a finer grid at points where the function changes more rapidly. Nov 8 '16 at 15:03
• @yohbs you do have two iterations of numerical approximation. But yes, Interpolation is the superior option in most cases, though I've encountered situations where it would fail. Nov 8 '16 at 15:08

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*)

• The robust solution is _?NumericQ, not _Real. None of 1, Pi, Sqrt[2], 1+I are Real, yet it's not 100% clear (though I admit likely) that NIntegrate could never plug these in. Nov 8 '16 at 15:26