Series expansion of InterpolatingFunction obtained from NDSolve

Question

I am trying to obtain a series expansion of the numerical solution of a differential equation. I encounter difficulties going beyond first-order expansions which I believe might be due to my inability to choose the right options in the functions involved

Consider as an easy example the following ODE for the exponential function:

Clear[y]
y[t_] = y[t] /.First@NDSolve[{y'[t] == y[t], y[0] == 1}, y[t], {t, 0, 1}]

We can check visually that the solution y[t] satisfies the equation y'[t]==y[t] and that it is equal to Exp[t] by following Wolfram's advice:

Plot[y'[t] - y[t] // RealExponent, {t, 0, 1}]

Plot[y[t] - Exp[t] // RealExponent, {t, 0, 1}]

This all seems very good.

I can then compute and compare the series expansions using

Series[Exp[t], {t, 0, 5}] // N // Chop//Normal
Series[y[t], {t, 0, 5}]//Normal

resulting in

(*
1. + 1. t + 0.5 t^2 + 0.166667 t^3 + 0.0416667 t^4 + 0.00833333 t^5
1. + 1. t + 2.00018 t^2 - 10892.9 t^3
*)

As can be seen, the first two coefficients coincide, and y[t] is only expanded up to order three.

From the documentation I learned that Interpolation fits polynomials between data points, which explains why y[t] only expands up to order three, and maybe -- depending on which data points are used to fit the polynomials -- the different coefficients as well.

To solve the problem, I tried setting the InterpolationOrder option in NDSolve to something larger than 3. This, however, resultet in Mma running out of memory and the kernel shutting down when trying to compute the series expansion of y[t].

I also read in the documentation that the option InterpolationOrder -> All 'specifies that the interpolation order should be chosen to be the same as the order of the underlying solution method', which suggests that the default underlying solution method may have an order which does not allow for an InterpolationOrder larger than 3.

Question: How can one obtain accurate series expansions of numerical solutions to differential equations up to arbitrary order?

Answer 1

The problem you're running into is that the usual InterpolatingFunction solutions from NDSolve have InterpolationOrder->3. This means that 4 derivatives of your result will be identically 0.

Now, it is possible to raise the order of the solutions by choosing different methods, but this approach is limited by what methods are available to be used. Instead, I would augment the NDSolve call so that it directly produces higher derivatives of your solution. For instance, given your initial NDSolve example:

NDSolveValue[{y'[t] == y[t], y[0] == 1}, y[x], {t, 0, 1}]

we could augment it by one order to get:

sol = NDSolveValue[
    {yp'[t] == yp[t], y'[t] == yp[t], y[0] == 1, yp[0] == 1},
    {y, yp},
    {t, 0, 2}
];

I used NDSolveValue because I find it more convenient. Also, note that I had to give the derivative a different name. Here's a comparison of the two solutions:

Plot[{sol[[1]]'[x] - Exp[x], sol[[2]][x] - Exp[x]}, {x, 0, 2}, PlotRange->All]

The interpolating function corresponding to yp gives a much more accurate result when compared to y'. However, augmenting the NDSolve/NDSolveValue call is a bit tedious, so here's a function that does this in general:

raiseOrder[Inactive[nds:NDSolve|NDSolveValue][eqns_, v_, {x_, r__}, opts:OptionsPattern[]], order_] := Module[
    {dpatt, eqn, ic, x0, max, low, neweqns},

    dpatt = Replace[v, List[d__]:>Alternatives[d]];

    (* Separate ODEs from initial/boundary conditions *)
    {eqn, ic} = Lookup[{False, True}] @ GroupBy[eqns, FreeQ[x]];

    (* Find initial/boundary points *)
    x0 = Cases[ic, (dpatt|Derivative[_][dpatt])[z_] :> z, Infinity] //Union;

    (* Boundary conditions are not supported *)
    If[MatchQ[x0, {_}],
        x0 = First @ x0,
        Return[Inactive[NDSolve][eqns, v, {x, r}, opts]]
    ];

    (* Find max order for each dependent variable v *)
    max = GroupBy[
        Cases[eqn, Derivative[n_][d:dpatt][x]:>{d, n}, Infinity],
        First->Last,
        Max
    ];

    (* Differentiate ODEs requested number of times *)
    neweqns = NestList[
        D[#, x]&,
        eqn,
        order
    ];

    (* Include needed additional first order equations. Also,
     * include substitutes for the lower order equations in neweqns *)
    low = KeyValueMap[
        Table[derivative[n-1][#1]'[x] == derivative[n][#1][x], {n, #2-1}]&,
        max + order
    ];

    (* Add initial conditions related to the additional ODEs in neweqns *)
    ic = Flatten @ Join[
        Most @ neweqns,
        ic
    ] /. Derivative->derivative /. x->x0;

    (* Evaluate NDSolve/NDSolveValue *)
    nds[
        Flatten @ Join[
            Take[neweqns, -1] /. Derivative[n_][y_]:>derivative[n-1][y]',
            low,
            ic
        ],
        Flatten @ KeyValueMap[
            Table[derivative[n][#1], {n, 0, #2+order-1}]&,
            max
        ],
        {x, r},
        opts
    ] /. derivative -> Derivative
]

derivative[0][y_] := y

Here is your example computed to order 10:

sol = First @ raiseOrder[
    Inactive[NDSolve][{y'[x]==y[x], y[0]==1}, y, {x, 0, 1}],
    10
]

Finally, here is the power series expanded around 0:

Series[y[x], {x, 0, 10}] /. sol //TeXForm

$1.`+x+0.5` x^2+0.16666666666666666` x^3+0.041666666666666664` x^4+0.008333333333333333` x^5+0.001388888888888889` x^6+0.0001984126984126984` x^7+0.0000248015873015873` x^8+\text{2.7557319223985893$\grave{ }$*${}^{\wedge}$-6} x^9+\text{2.755731922398589$\grave{ }$*${}^{\wedge}$-7} x^{10}+O\left(x^{11}\right)$

Compare this to the actual series:

Series[Exp[1. x], {x, 0, 10}] //TeXForm

$1+x+0.5` x^2+0.16666666666666666` x^3+0.041666666666666664` x^4+0.008333333333333333` x^5+0.0013888888888888887` x^6+0.00019841269841269839` x^7+0.000024801587301587298` x^8+\text{2.7557319223985884$\grave{ }$*${}^{\wedge}$-6} x^9+\text{2.7557319223985883$\grave{ }$*${}^{\wedge}$-7} x^{10}+O\left(x^{11}\right)$

Of course, this isn't a very good example because the expansion around 0 doesn't really require solving an ODE. Instead, let's compute the power series around 1:

Series[y[x], {x, 1, 10}] /. sol //TeXForm

$2.718281861392969`+2.718281861392969` (x-1)+1.3591409306964846` (x-1)^2+0.45304697689882817` (x-1)^3+0.11326174422470704` (x-1)^4+0.022652348844941408` (x-1)^5+0.0037753914741569016` (x-1)^6+0.0005393416391652716` (x-1)^7+0.00006741770489565895` (x-1)^8+\text{7.490856099517662$\grave{ }$*${}^{\wedge}$-6} (x-1)^9+\text{7.490856099517661$\grave{ }$*${}^{\wedge}$-7} (x-1)^{10}+O\left((x-1)^{11}\right)$

And the actual answer:

Series[Exp[x], {x, 1., 10}] //TeXForm

$2.718281828459045`+2.718281828459045` (x-1.`)+1.3591409142295225` (x-1.`)^2+0.45304697140984085` (x-1.`)^3+0.11326174285246021` (x-1.`)^4+0.02265234857049204` (x-1.`)^5+0.0037753914284153404` (x-1.`)^6+0.0005393416326307629` (x-1.`)^7+0.00006741770407884536` (x-1.`)^8+\text{7.490856008760597$\grave{ }$*${}^{\wedge}$-6} (x-1.`)^9+\text{7.490856008760595$\grave{ }$*${}^{\wedge}$-7} (x-1.`)^{10}+O\left((x-1.`)^{11}\right)$

Answer 2

Using Method->"ExplicitRungeKutta" with a larger value of the option "DifferenceOrder" allows recovering more terms of the series expansion.