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

Question

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?

Answer 1

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 DifferentialEquations`InterpolatingFunctionAnatomy` 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}]]}]

Answer 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}]]}]]

Answer 3

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}]]}]