7
$\begingroup$

For example given the code:

s = NDSolve[{y''[x] + Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0}, 
  y, {x, 0, 30}]
Plot[Evaluate[{y[x], y'[x], y''[x]} /. s], {x, 0, 30}, 
 PlotStyle -> Automatic]

giving output enter image description here

Are the derivatives shown here the numerical derivatives the NDSolve would need to calculate when solving the differential equation, or are the just derivatives determined after the fact from the spline of the interpolating function?

My expectation is that these are just derivatives of the spline but it would be nice if these were the values calculated by NDSolve.

Improve this question
$\endgroup$
2
  • $\begingroup$ "it would be nice if these were the values calculated by NDSolve." - then you need to actually provide the relevant equations for the derivatives, since y'[x] is by necessity directly differentiating the InterpolatingFunction. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 23, 2022 at 15:25
  • 1
    $\begingroup$ Interpolation can construct an InterpolatingFunction that reproduces the values of derivatives as well as those of the function itself. Since this capability exists, I would think that NDSolve would provide that information in its InterpolatingFunction output. But this is speculation on my part. As an aside, you could use NDSolveValue[yourEquation, {y, y', y''}, {x, 0, 30}] to get separate InterpolatingFunction objects representing the function and each of its derivatives, right from within NDSolve. $\endgroup$
    – MarcoB
    Commented May 23, 2022 at 15:27

1 Answer 1

Reset to default
7
$\begingroup$

The default type of interpolation produced by NDSolve stores for each discrete step the domain abscissa $x_j$, the function value $y(x_j)$, and the derivative values $y^{k}(x_j)$ up to the order of the ODE. In the OP's example, the derivatives stored are $y'(x_j)$ and $y''(x_j)$. The data structure inside InterpolatingFunction[<...>] is shown in What's inside InterpolatingFunction[{{1., 4.}}, <>]? and in the code toward the end of this answer, Interpolating data with a step. Part 4 contains the function and derivative data for each step. It is stored in a flat array of the form

(*
{y[x0], y'[x0],..., Derivative[k][y][x0],
 y[x1], y'[x1],..., Derivative[k][y][x1],
 ...}
*)

The data also contains offsets that indicate where each function/derivative value unit ends. When the order of the ODE is k = 2, there are three values per unit y, y', y''.

yIFN = y /. First[s];
yIFN[[4]]
(*
{Developer`PackedArrayForm,       (* type of data *)
 {0, 3, 6, 9, ... 1407},          (* offsets *)
 {1.`, 0.`, -0.8414709848078965`, (* function/derivative values *)
  0.9999999750347968`, -0.00014493961665665655`, -0.8414709503116456`,
  0.9999999251043933`, -0.00028987922142966306`, -0.8414708813191469`,
  ...,
  0.9989990862203024`,  0.04101857645543858`,    -0.8400880807981383`}}
*)

The argument "ValuesOnGrid" gives a slice of this data for each of the derivatives computed by NDSolve (i.e., y, y', y'' in the OP example):

yIFN["ValuesOnGrid"][[{1, 2, 3, -1}]]
yIFN'["ValuesOnGrid"][[{1, 2, 3, -1}]]
yIFN''["ValuesOnGrid"][[{1, 2, 3, -1}]]

(*       y[x]         |           y'[x]            |         y''[x]
{1.`,                 | {0.`,                      | {-0.8414709848078965`,
 0.9999999750347968`, |  -0.00014493961665665655`, |  -0.8414709503116456`,
 0.9999999251043933`, |  -0.00028987922142966306`, |  -0.8414708813191469`, 
 0.9989990862203024`} |  0.04101857645543858`}     |  -0.8400880807981383`} 
*)

Note these values match the data shown for yIFN[[4]]. If you ask for higher derivatives, they are computed from the interpolation. We have 6 values total for each interpolation interval (3 at each end point). They determine a degree 5 interpolant. Derivatives of order 6 and higher are zero.

Derivative[3][yIFN]["ValuesOnGrid"][[{1, 2, 3, -1}]]
Derivative[4][yIFN]["ValuesOnGrid"][[{1, 2, 3, -1}]]
Derivative[5][yIFN]["ValuesOnGrid"][[{1, 2, 3, -1}]]
Derivative[6][yIFN]["ValuesOnGrid"][[{1, 2, 3, -1}]]
(* 
{-146559., -146559., -146559., ..., -0.0565191}
{5.10523*10^9, -5.10523*10^9, -5.10523*10^9, ..., 1.1701}
{-5.92786*10^13, -5.92786*10^13, -5.92786*10^13, ..., 0.366953}
{0., 0., 0., ..., 0.}
*)

What looks like numerical instability is probably due to the small step sizes at the beginning:

Differences[yIFN["Grid"]][[{1, 2, 3, -1}]]
(*  {{0.000172246}, {0.000172246}, {0.00382697}, {0.0971486}}  *)

Derivative[5][yIFN] // ListLinePlot[#, PlotRange -> All] &
Derivative[5][yIFN] // ListLinePlot

Update

Another thing worth mentioning is that only one thing is changed in yIFN when it is differentiated. The order of the derivative is changed in the InterpolatingFunction structure.

yIFN''''[[4]] === yIFN[[4]]
(*  True  *)

The derivative of the interpolation may be inspected by Part or via the "DerivativeOrder" method:

yIFN[[2, 6]]
yIFN''[[2, 6]]
yIFN''''[[2, 6]]
(*
  0
  {2}
  {4}
*)

yIFN''''["DerivativeOrder"]
(*  Derivative[4]  *)

Related Q&A:

These deal with variations due to InterpolationOrder -> All:

Improve this answer
$\endgroup$
6
  • $\begingroup$ Does changing InterpolationOrder when calling NDSolve cause Mathematica to store higher-order derivatives in the InterpolatingFunction? $\endgroup$
    – Michael Seifert
    Commented May 23, 2022 at 15:36
  • 3
    $\begingroup$ @MichaelSeifert InterpolationOrder -> All will change things; but a lower InterpolationOrder ->3 seems to have no effect; and a higher InterpolationOrder -> 7 does the same thing as All. InterpolationOrder -> All will store an interpolant for each step of a degree equal to the order of the NDSolve method. For the default LSODA, which has a variable order, it uses either Hermite order-3 or local Taylor series. For Runge-Kutta methods, it uses Chebyshev series, which can only be used (in InterpolatingFunction) if the order is fixed. I left these complications out of my answer. $\endgroup$
    – Michael E2
    Commented May 23, 2022 at 15:48
  • $\begingroup$ IIRC Method -> "Extrapolation" also forces NDSolve[] to use a Chebyshev interpolant behind the scenes. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 23, 2022 at 15:54
  • 1
    $\begingroup$ @J.M. The method(s) used (for each step) may be inspected with yIFN[[5]], but I didn't want to go through all the methods. $\endgroup$
    – Michael E2
    Commented May 23, 2022 at 15:59
  • $\begingroup$ Thank you for the very detailed answer @MichaelE2! This helps a ton! $\endgroup$
    – akozi
    Commented May 23, 2022 at 17:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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