3
$\begingroup$

Start with

H = q[t]^2 + p[t]^2 + q[t]^3
A = D[H, p[t]]
B = -D[H, q[t]]
NDSolve[{q'[t] == A, p'[t] == B, q[0] == p[0] == .2}, {q[t], 
    p[t]}, {t, 0, 5}] // First;
R = Evaluate[{q[t], p[t]} /. %]

This gives

>> {InterpolatingFunction[{{0., 5.}}, <>][t], InterpolatingFunction[{{0., 5.}}, <>][t]}

Which looks like

ParametricPlot[R, {t, 0, 4}, PlotStyle -> Red]

enter image description here.

Now, take 

Series[R, {t, 0, 10}] // Normal

which gives

{0.2 + 0.4 t - 2.08063 t^2 + 8477.23 t^3, 
 0.2 - 0.52 t - 2.5593 t^2 + 10427.5 t^3} .

Something is utterly wrong, here. This should be an infinite polynomial. I don't even have to plot it to know it won't, obviously, approximate my original interpolating function correctly.

Edit:

  1. Is there a maximum number of derivatives one can take from an interpolating function? Why the sudden truncation? - Answered in the comments.
  2. Is there a way to circumvent this?
|improve this question
$\endgroup$
  • 7
    $\begingroup$ An InterpolatingFunction interpolates between points precicely by using interpolating polynomials. The default order seems to be 3 (but I don't know if there are cases where Mathematica automatically uses a higher or lower order), explaining why your Series truncate after t^3. $\endgroup$ – Marius Ladegård Meyer Dec 14 '16 at 20:43
  • 1
    $\begingroup$ For an answer to your second question: probably not. $\endgroup$ – march Dec 14 '16 at 20:47
  • 1
    $\begingroup$ Using InterpolationOrder -> All in NDSolve will yield an interpolation with an order equal to the order of the step (which is variable with the default method). It still won't be of order 10. $\endgroup$ – Michael E2 Dec 14 '16 at 23:38
  • 1
    $\begingroup$ Using the Bulirsch-Stoer method (Method -> "Extrapolation") in conjunction with Michael's advice should net you a few more coefficients, but maybe not as many as you'd like. $\endgroup$ – J. M. will be back soon♦ Dec 15 '16 at 1:42
  • 1
    $\begingroup$ @QuantumBrick InterpolationOrder is in the base documentation of NDSolve[] and such this question is at risk of closure. There have also been multiple votes not to close this question as the solution isn't as obvious as many others that get closed for this reason. $\endgroup$ – Feyre Dec 15 '16 at 17:57
6
$\begingroup$

You might use FunctionInterpolation to re-interpolate to a higher order: 

g = FunctionInterpolation[R[[1]], {t, 0, 4}, InterpolationOrder -> 6]

then you get a higher order series approximation:

s[t_] = Series[g[t], {t, 0, 6}] // Normal

0.2 + 0.456923 t - 0.780614 t^2 - 0.132724 t^3 + 0.325397 t^4 - 0.078419 t^5

Show[{
  Plot[R[[1]], {t, 0, 4}, PlotStyle -> Red], Plot[s[t], {t, 0, 4}]}, 
   PlotRange -> {-1.5, 1}]

enter image description here

note the cubic series (green line) is really poor, since its only good over the first increment from NDSolve which is ~0.00012

the parametric form:

s[t_] = {Series[FunctionInterpolation[#, {t, 0, 4},
        InterpolationOrder -> 10][t], {t, 0, 20}] // Normal} & /@ R
Show[{
  ParametricPlot[R, {t, 0, 3.546}, PlotStyle -> {Thickness[.02], Red}],
  ParametricPlot[s[t], {t, 0, 3.546}, PlotStyle -> Blue]}]

enter image description here

|improve this answer
$\endgroup$
3
$\begingroup$

There are higher order methods available in NDSolve, for example "ImplicitRungeKutta".

H = q[t]^2 + p[t]^2 + q[t]^3;
A = D[H, p[t]];
B = -D[H, q[t]];
NDSolve[{q'[t] == A, p'[t] == B, q[0] == p[0] == .2},
   {q[t], p[t]}, {t, 0, 5}, 
   Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}, (* for degree 10 polys *) 
   InterpolationOrder -> All] // First;
R = Evaluate[{q[t], p[t]} /. %];

Series[R, {t, 0, 10}] // Normal
(*
  {0.2 + 0.4 t - 0.52 t^2 - 0.425533 t^3 - 97.2464 t^4 + 
    3.83907*10^6 t^5 - 8.07936*10^10 t^6 + 9.74544*10^14 t^7 - 
    6.7802*10^18 t^8 + 2.5365*10^22 t^9 - 3.95812*10^25 t^10, 
   0.2 - 0.52 t - 0.64 t^2 + 0.392703 t^3 + 169.596 t^4 - 
    6.65809*10^6 t^5 + 1.40133*10^11 t^6 - 1.69034*10^15 t^7 + 
    1.17602*10^19 t^8 - 4.39952*10^22 t^9 + 6.86524*10^25 t^10}
*)

This has significant error in the higher degree terms. In part at least, this is due to expanding at t == 0, where NDSolve takes a few very small steps as it gets its bearings. It behaves better in the interior:

ser = Series[R, {t, 1, 10}] // Normal
(*
  {-0.0105103 - 0.592927 (-1 + t) + 0.0206892 (-1 + t)^2 + 
     0.382821 (-1 + t)^3 - 0.18246 (-1 + t)^4 - 0.0667898 (-1 + t)^5 + 
     0.114268 (-1 + t)^6 - 0.0269471 (-1 + t)^7 - 0.0295921 (-1 + t)^8 - 
     0.0847381 (-1 + t)^9 + 1.31425 (-1 + t)^10,
   -0.296464 + 0.0206892 (-1 + t) + 0.574232 (-1 + t)^2 - 0.364921 (-1 + t)^3 - 
     0.166974 (-1 + t)^4 + 0.342808 (-1 + t)^5 - 0.0945715 (-1 + t)^6 - 
     0.123299 (-1 + t)^7 + 0.129194 (-1 + t)^8 - 2.51112 (-1 + t)^9 + 36.9119 (-1 + t)^10}
*)

The polynomials are accurate only over a few steps near the center of the series. The method Method -> "Extrapolation" does quite a bit better, but sometimes with a lower degree than 10 (such as 8 or 9); however, it's still not accurate over the whole curve as george2079's method seems to be. But one should keep in mind that the usual purpose of series is local approximation.

|improve this answer
$\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.