# How can I know how many points are used to by the interpolating function returned by NDSolve?

In a problem in which the solution has oscillations that I solved, the amplitude of the oscillations is much lower than I expected. It looks as if the interpolating function doesn't use enough points to plot the results correctly. Has anyone encountered a problem like this before?

It would help a lot if someone could show me a way to deal with this.

The left image shows what I expected, and the right one is what I get.

The equations are from this paper (also here).

The following is my code; it is a little bit long...

f[x_] = x;
n = 50;
\[CapitalEpsilon]\[CapitalEpsilon] = {0, 5.2, 12.4, 19.1, 23.8, 31.3,
50};
\[CapitalEpsilon] =
DeleteDuplicates[
Flatten[Table[
N[Rationalize[
Array[ f,
Round[(n (\[CapitalEpsilon]\[CapitalEpsilon][[i +
1]] - \[CapitalEpsilon]\[CapitalEpsilon][[i]]))/
50], {\[CapitalEpsilon]\[CapitalEpsilon][[i]], \
\[CapitalEpsilon]\[CapitalEpsilon][[i + 1]]}]]], {i, 1, 6}]]];
n = Length[\[CapitalEpsilon]];
\[ScriptCapitalD] = 1/2 (1 - Sqrt[1 - (10/t)^2])^2;

For[j = 2, j <= n, j++,
P[\[CapitalEpsilon][[j]], t] = {Px[\[CapitalEpsilon][[j]], t],
Py[\[CapitalEpsilon][[j]], t], Pz[\[CapitalEpsilon][[j]], t]}]
For[j = 2, j <= n, j++,
Pbar[\[CapitalEpsilon][[j]], t] = {Pbarx[\[CapitalEpsilon][[j]], t],
Pbary[\[CapitalEpsilon][[j]], t], Pbarz[\[CapitalEpsilon][[j]], t]}]

B = {0.02, 0, 0.9998};

\[Zeta] = 1.202;
Ebare = 10;
Ebaree = 15;
Ebarx = 24;
\[Beta]e = 0.315;
\[Beta]ee = 0.21;
\[Beta]x = 0.131;
u = 10^51*(0.197*10^-18)^2*1000*Sqrt[2]*1.16638*10^-5/(3*10^5);
w = 1/0.197;

PP = Table[
D[P[\[CapitalEpsilon][[i]], t], t] ==
w/\[CapitalEpsilon][[i]] Cross[B, P[\[CapitalEpsilon][[i]], t]] +
u \[ScriptCapitalD]*
Cross[Sum[
2/(3 \[Zeta]) (((\[Beta]e (\[Beta]e*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]e*\[CapitalEpsilon][[i]]) +
1) Ebare) + (\[Beta]x (\[Beta]x*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]x*\[CapitalEpsilon][[i]]) +
1) Ebarx)) P[\[CapitalEpsilon][[i]],
t] - ((\[Beta]ee (\[Beta]ee*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]ee*\[CapitalEpsilon][[i]]) +
1) Ebaree) + (\[Beta]x \
(\[Beta]x*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]x*\[CapitalEpsilon][[i]]) +
1) Ebarx)) Pbar[\[CapitalEpsilon][[i]],
t]) (\[CapitalEpsilon][[i]] - \[CapitalEpsilon][[i -
1]]), {i, 2, n}], P[\[CapitalEpsilon][[i]], t]], {i, 2,
n}];
PV = Table[
P[\[CapitalEpsilon][[i]], t] == {0,
0, (-((0.00009367045833333334 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.131 \[CapitalEpsilon][[i]]))) + (
0.0031255875000000006 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.315 \[CapitalEpsilon][[i]])))/((
0.00009367045833333334 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.131 \[CapitalEpsilon][[i]])) + (
0.0031255875000000006 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.315 \[CapitalEpsilon][[i]])))} /. t -> 10, {i, 2, n}];
PbarP = Table[
D[Pbar[\[CapitalEpsilon][[i]], t],
t] == -w/\[CapitalEpsilon][[i]] Cross[B,
Pbar[\[CapitalEpsilon][[i]], t]] +
u \[ScriptCapitalD]*
Cross[Sum[
2/(3 \[Zeta]) (((\[Beta]e (\[Beta]e*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]e*\[CapitalEpsilon][[i]]) +
1) Ebare) + (\[Beta]x (\[Beta]x*\[CapitalEpsilon][[i]])^2)/( (E^(\[Beta]x*\[CapitalEpsilon][[i]]) +
1) Ebarx)) P[\[CapitalEpsilon][[i]],
t] - ((\[Beta]ee (\[Beta]ee*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]ee*\[CapitalEpsilon][[i]]) +
1) Ebaree) + (\[Beta]x \
(\[Beta]x*\[CapitalEpsilon][[
i]])^2)/( (E^(\[Beta]x*\[CapitalEpsilon][[i]]) +
1) Ebarx)) Pbar[\[CapitalEpsilon][[i]],
t]) (\[CapitalEpsilon][[i]] - \[CapitalEpsilon][[i -
1]]), {i, 2, n}], Pbar[\[CapitalEpsilon][[i]], t]], {i,
2, n}];
PbarV = Table[
Pbar[\[CapitalEpsilon][[i]], t] == {0,
0, (-((0.00009367045833333334 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.131 \[CapitalEpsilon][[i]]))) + (
0.0006173999999999999 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.21 \[CapitalEpsilon][[i]])))/((
0.00009367045833333334 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.131 \[CapitalEpsilon][[i]])) + (
0.0006173999999999999 \[CapitalEpsilon][[i]]^2)/(
1 + E^(0.21 \[CapitalEpsilon][[i]])))} /. t -> 10, {i, 2,
n}];
s = Join[PP, PbarP, PV, PbarV];

sss = NDSolve[s,(*Z*)Pz[\[CapitalEpsilon][[5]], t], {t, 10, 200}]
Plot[Pz[\[CapitalEpsilon][[5]], t] /. sss, {t, 10, 200},
AxesOrigin -> {10, 0}, PlotRange -> {{10, 200}, {-1, 1}},
Frame -> True]


If anyone want to know more clearly about my question(what equation I am solving) and find the paper hard to read, I will show some other code below.

  P = {Px[\[CapitalEpsilon], t], Py[\[CapitalEpsilon], t],
Pz[\[CapitalEpsilon], t]};
Pbar = {Pbarx[\[CapitalEpsilon], t], Pbary[\[CapitalEpsilon], t],
Pbarz[\[CapitalEpsilon], t]};

D[P, t] = {+\[Omega][\[CapitalEpsilon]] B + \[Sqrt]2 Gf \
Integrate[((ne[\[CapitalEpsilon], r] +
nx[\[CapitalEpsilon], r]) P - (nebar[\[CapitalEpsilon],
r] + nx[\[CapitalEpsilon],
r]) Pbar), \[CapitalEpsilon]]}\[Cross]P;

D[Pbar, t] = {-\[Omega][\[CapitalEpsilon]] B + \[Sqrt]2 Gf \
Integrate[((ne[\[CapitalEpsilon], r] +
nx[\[CapitalEpsilon], r]) P - (nebar[\[CapitalEpsilon],
r] + nx[\[CapitalEpsilon],
r]) Pbar), \[CapitalEpsilon]]}\[Cross]Pbar;


\[Omega], B, Gf, ne[\[CapitalEpsilon], r], nx[\[CapitalEpsilon], r], \ nebar[\[CapitalEpsilon], r], nxbar[\[CapitalEpsilon], r] are known.

In my origin code before, I try to use discrete \[CapitalEpsilon] instead of continuous \[CapitalEpsilon]

• 1. This question is unanswerable without seeing your code. 2. Try increasing the PlotPoints setting in Plot[], as a first thing to try. Oct 10, 2018 at 4:49
• Ok, I will show my whole code, but I think the question is the interpolatingfunction solved by NDSolve doesn't have enough points in oscillation.
– 袁子奕
Oct 10, 2018 at 7:34
• The For[] bothers me a bit, and I suspect this can be reformulated better. Where did these equations come from? Do you have a book or paper as a source for this? Oct 10, 2018 at 8:10
• @J.M.issomewhatokay. I'm a rookie of MMA, the equations come from a paper of Supernova neutrino, I try to repeat their result as shown in the left picture I upload, but finally what I get is the right one. Because they looks same, so I doubt whether the NDSolve give enough points in the places where oscillate rapidly.
– 袁子奕
Oct 10, 2018 at 8:38
• ……The paper is 27 pages long and there're 85 equations listed in that paper, but only a few of them are solved in your code, right? Then, which ones? You mentioned you've discretized the original equation all by yourself, then, in what way? These should be clarified in your question, or we can't even check if the equation is "translated" correctly. Oct 10, 2018 at 13:27

Interpolation[...][[-3,1]]


returns the time values of the interpolation-object.

sol = NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t]^3, x[0] == y[0] == 1},
{x, y}, {t, 20}]


{{x->InterpolatingFunction[{{0.,20.}},<>],y->InterpolatingFunction[{{0.,20.}},<>]}}

Length[(x /. sol[[1]])@"Grid"]


298

Also

Length[(x /. sol[[1]])@"ValuesOnGrid"]


298

Both match the result from the method suggested in Ulrich Neumann's answer:

Length[(x /. sol[[1]])[[3, 1]]]


298

• Interesting answer, I didn't know the option "Grid", where/how did you find it? Oct 10, 2018 at 7:44
• @UlrichNeumann, I have seen "ValuesOnGrid" somewhere on this site some time back. "Grid" was lucky guess.
– kglr
Oct 10, 2018 at 7:49
• do you think 'StartingStepSize -> 1/10, Method -> {"FixedStep", Method -> "ExplicitEuler"}' will help to make the InterpolatingFunction have more points?
– 袁子奕
Oct 10, 2018 at 8:02
• @袁子奕, I don;t know how to control the number of points sampled by NDSolve. I would first try a large value for the option PlotPoints in Plot as suggested by J. M.issomewhatokay.
– kglr
Oct 10, 2018 at 8:07
• @袁子奕 Notice the initial value problem (IVP) solver of NDSolve is very robust, so does Plot. So in most cases an unexpected result only suggests one thing: the equations system itself is wrong. Oct 10, 2018 at 13:31