5
$\begingroup$
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?

$\endgroup$
6
$\begingroup$

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

$\endgroup$
4
$\begingroup$

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

$\endgroup$
  • $\begingroup$ It works well, thanks! $\endgroup$ – SunnySky Sep 18 '12 at 8:03
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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