Sudden truncation of numerical series expansion

Question

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]

.

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?

Answer 1

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

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

Answer 2

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.