4
$\begingroup$

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

enter image description here

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

enter image description here

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ NDSolve[{y'[t] == y[t], y[0] == 1}, y'[t], {t, 0, 1}] will solve for y'[t]. By solving for arbitrary derivatives, you might be able to use Taylor series expansion directly? $\endgroup$
    – user1722
    Oct 30, 2014 at 13:37
  • $\begingroup$ Thanks @barrycarter. But y''[t] /. First@NDSolve[{y'[t] == y[t], y[0] == 1}, y''[t], {t, 0, 1}] /. t -> 0 returns 4.00037, not one, as I would like it to... $\endgroup$
    – Eckhard
    Oct 30, 2014 at 13:43
  • $\begingroup$ I believe trying to get orders higher than the ODE itself is doomed. Can't you try to fit the points with a nice model? $\endgroup$
    – Dr. belisarius
    Oct 30, 2014 at 13:46
  • $\begingroup$ Do y''[t] /. First@NDSolve[{y'[t] == y[t], y[0] == 1}, y''[t], {t, 0, 1}] and then Plot[%,{t,0,1}] and the result looks correctish (I think your t->0 is being evaluated at the wrong place)? I agree with @Eckhard too: try FindSequenceFunction or something maybe. OK, I see what happens. The function is 4+ at t=0 (PlotRange->All helped) but it becomes correct even an infinitesimal distance from t=0. That's sort of reasonable since t=0 is a boundary. $\endgroup$
    – user1722
    Oct 30, 2014 at 14:26
  • $\begingroup$ Using Method->"ExplicitRungeKutta" with a large enough "DifferenceOrder" seems to do the trick. I might have to brush up on my numerical analysis skills to be better able to predict the effect of various methods and options. $\endgroup$
    – Eckhard
    Oct 30, 2014 at 15:10

2 Answers 2

Reset to default
4
$\begingroup$

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]

enter image description here

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
]

enter image description here

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)$

Improve this answer
$\endgroup$
3
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Method -> {"FixedStep", Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 5}} is also available, there should be more, I guess. $\endgroup$
    – xzczd
    Nov 13, 2014 at 4:35
  • $\begingroup$ Note that for "ImplicitRungeKutta", the difference order can be raised to arbitrarily large values (subject to the normal constraints of numerical precision, memory, etc.). The default integration is based on Gaussian integration, and the solutions should converge rapidly as a function of order for ODEs with analytic coefficients. The coefficients up to about order order/2 for the example in the OP. Might be good to reiterate that InterpolationOrder -> All is needed to get all the terms. $\endgroup$
    – Michael E2
    Nov 16, 2017 at 13:51

Your Answer

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

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