# How could I get the value of y[t] at each specific interpolation point?

sol = NDSolve[{Derivative[2][y][t] + Sin[y[t]] == 0, Derivative[1][y][0] == 0, y[0] == 1},
y, {t, 0, 2}]


the above-mentioned differential equations can be solved by Numerical interpolated integrations at a series of discontinuous points, and the number of the points depends on the Integration methods selected by the NDSolve. My problem is that how I can obtain all these points selected because I want to get the specific value of y[t] at each interpolation point?

-

To see the points being kept around by the InterpolatingFunction[] object spat out by NDSolve[], one could use the special syntax for extracting the "anatomy" of an InterpolatingFunction[] object (used internally by the functions in the DifferentialEquationsInterpolatingFunctionAnatomy package).

As a concrete example:

sol = y /. First@NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 1, y'[0] == 0}, y, {t, 0, 2}];

xa = First @ sol["Coordinates"]
{0., 0.00010487, 0.000209741, 0.00487198, 0.00953422, 0.0141965,
0.0142281, 0.0142597, 0.0142913, 0.0143229, 0.0143545, 0.0144177,
0.014481, 0.0145442, 0.0151764, 0.0158087, 0.0164409, 0.0170732,
0.0233957, 0.0297182, 0.0360407, 0.0423632, 0.0486857, 0.111911,
0.175136, 0.238361, 0.301586, 0.364811, 0.428036, 0.498301, 0.568567,
0.638832, 0.709098, 0.779363, 0.849629, 0.919894, 1.0109, 1.10191,
1.19292, 1.28393, 1.37494, 1.44965, 1.52436, 1.59907, 1.67378,
1.74849, 1.8232, 1.89792, 1.94896, 2.}

ya = sol["ValuesOnGrid"]
{1., 1., 1., 0.99999, 0.999962, 0.999915, 0.999915, 0.999914,
0.999914, 0.999914, 0.999913, 0.999913, 0.999912, 0.999911, 0.999903,
0.999895, 0.999886, 0.999877, 0.99977, 0.999628, 0.999454, 0.999245,
0.999003, 0.994734, 0.987113, 0.976157, 0.961891, 0.944346, 0.923564,
0.896729, 0.866039, 0.831588, 0.793491, 0.751875, 0.706889, 0.658698,
0.591821, 0.520327, 0.44472, 0.365554, 0.283431, 0.214255, 0.143892,
0.0727304, 0.00116316, -0.0704105, -0.141592, -0.211986, -0.259428,
-0.306201}

Plot[sol[x], {x, 0, 2}, Epilog -> {AbsolutePointSize[5], Red, Point[Transpose[{xa, ya}]]}]


-

Alternative way (not general, would have to be tweaked if you change InterpolationOrder ...) :

x-values used to define the interpolation :

 xdata = sol[[1, 1, 2]][[3, 1]];


{y,y',y''}-values :

ydataRaw = sol[[1, 1, 2]][[4, 3]];

ydata = ydataRaw[[Range[1, Length[ydataRaw], 3]]];
ypdata = ydataRaw[[Range[2, Length[ydataRaw], 3]]];
yppdata = ydataRaw[[Range[3, Length[ydataRaw], 3]]];


data to be plotted :

data = Transpose[{xdata, ydata}];
datap = Transpose[{xdata, ypdata}];
datapp = Transpose[{xdata, yppdata}];

GraphicsColumn[{Show[ListPlot[data, PlotStyle -> Red], Plot[y[t] /. sol[[1]], {t, 0, 2}]],
Show[ListPlot[datap, PlotStyle -> Red], Plot[y'[t] /. sol[[1]], {t, 0, 2}]],
Show[ListPlot[datapp, PlotStyle -> Red], Plot[y''[t] /. sol[[1]], {t, 0, 2}]]}]


-

Using the examples from the docs on EvaluationMonitor

  {sol, evaluations} =  Reap[NDSolve[{Derivative[2][y][t] + Sin[y[t]] == 0,
Derivative[1][y][0] == 0, y[0] == 1}, y, {t, 0, 2},
EvaluationMonitor :> Sow[{t, y[t], y'[t]}]]];

Show[Plot[Evaluate[First[{y[t], y'[t]} /. sol]], {t, 0, 5}],
ListPlot[{evaluations[[1, All, {1, 2}]], evaluations[[1, All, {1, 3}]]}]]


 evaluations


(Torn image made with Heike's torn function)

With a different setting for Method:

  {sol, evaluations} =  Reap[NDSolve[{Derivative[2][y][t] + Sin[y[t]] == 0,
Derivative[1][y][0] == 0, y[0] == 1}, y, {t, 0, 2},
EvaluationMonitor :> Sow[{t, y[t], y'[t]}],
Method -> "Extrapolation"]];

Show[Plot[Evaluate[First[{y[t], y'[t]} /. sol]], {t, 0, 5}],
ListPlot[{evaluations[[1, All, {1, 2}]], evaluations[[1, All, {1, 3}]]}]]


-
It works well, thanks! – SunnySky Sep 18 '12 at 8:03