Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
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?

share|improve this question

3 Answers 3

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

plot of approximate solution and interpolation points

share|improve this answer

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

enter image description here

 evaluations

enter image description here

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

enter image description here

share|improve this answer
    
It works well, thanks! –  SunnySky Sep 18 '12 at 8:03

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

plot

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.