Evaluate function defined by DifferentialRoot

Question

I have the following sequence of rationals that I want to find the generating function of:

d1 = {1, 1, 0, -(1/3), -(1/12), 3/20, 7/60, -(1/35), -(271/3360), -(31/1440), 2953/75600, 54941/1663200, -(32377/3991680), -(313627/12355200), -(1562149/201801600), 57637519/4540536000, 81128389/6918912000, -(5841262799/2470051584000), -(16908002893/1852538688000), -(122799753991/38398074624000), 38177259047/8439137280000, 73864046480467/16127191342080000, -(422047957817/622452999168000), -(4070047619579/1156513438080000), -(899675774837639/640028142673920000), 9035541454237687/5333567855616000000, 21789055343201563/11556063687168000000, -(139475898744014537/802320993137664000000), -(5703599712765137243/4032177298845696000000), -(100308434035944184283/157254914654982144000000), 580014783870860474531/888388154219704320000000}

FindSequenceFunction gives a DifferenceRoot object that can be evaluated to recover the values of the sequence.

sf = FindSequenceFunction[d1]
(*DifferenceRoot[Function[{y, n}, {(-2 - n^2) y[n] + 3 (1 + n + n^2) y][1 + n]] - 3 (1 + n)^2 y][2 + n]] +  3 (1 + n) (2 + n]) y][3 + n]] == 0, y[1] == 1, y][2] == 1, y][3] == 0}]]*)

By the way, why does Mathematice use these \[FormalN] variables in output like this?

sf[5]==-1/12
(*True*)

Similarly, FindGeneratingFunction returns a DifferentialRoot object:

gf[x_] = FindGeneratingFunction[d1, x]
(*DifferentialRoot[Function[{y,x},{(-6+3 x-3 x^2+2 x^3) y[x]+x (6-3 x+x^3) (y'[x]+x^2 (-3+3 x-3 x^2+x^3) y''[x]==0,y[0]==0,y'[0]==1}]][x]/x*)

However, evaluating gf does not seem to work. Also, when I try to integrate the differential equation found by FindGeneratingFunction directly using NDSolve the input is returned unevaluated.

I can obtain an approximation to the generating function as

f1[x_] = Dot[d1, x^Range[0, Length[d1] - 1]]

but this is of course only valid for small x and I'd prefer having a differential equation that I can solve numerically to working with truncated power series of large degree.

Question: How can I evaluate the DifferentialRoot object returned by FindGeneratingFunction in the example described above?

Answer 1

However, evaluating gf does not seem to work

When Mathematica returns DifferentialRoot as solution to a differential equation, it really means M could not solve the ODE. This is not the necessarily the case with DifferenceRoot as you found out (I do not know much about DifferenceRoot, but I've seen DifferentialRoot before many times). Here is an example:

ode = -y[x] + (1 + x)*Derivative[2][y][x] + x^2*Derivative[3][y][x] == 0;
DSolve[ode, y[x], x]

You can't use y[x] generated above. The above really means M could not solve the ODE. For your case, the ODE is

ode = (-6 + 3*x - 3*x^2 + 2*x^3)*y[x] + x*(6 - 3*x + x^3)*Derivative[1][y][x] + 
          x^2*(-3 + 3*x - 3*x^2 + x^3)*Derivative[2][y][x] == 0;
ic = {y[0] == 0, y'[0] == 1};
DSolve[{ode, ic}, y[x], x];

However, if you want the solution to the above ODE, Maple 18.1 gives this very complicated solution in terms of HeunG function.

restart;
ode := (-6 + 3*x - 3*x^2 + 2*x^3)*y(x) + x*(6 - 3*x + x^3)*diff(y(x),x) + 
         x^2*(-3 + 3*x - 3*x^2 + x^3)*diff(y(x),x$2)= 0;
ic := y(0)= 0, D(y)(0)= 1:
sol:=dsolve({ode, ic},y(x));

gives

y(x) = (-HeunG((2*I)*sqrt(3)/(I*sqrt(3)-3), 864/((I*sqrt(3)+3)^3*
    (I*sqrt(3)-3)^2), 1+I*sqrt(2), (72*I)*(sqrt(6)-(3*I)*sqrt(2)-(2*I)*sqrt(3)-
    6)*sqrt(2)/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1, 1, (-1+I*sqrt(3)+2^(2/3))
    /(I*sqrt(3)-3))*_C2/HeunG((2*I)*sqrt(3)/(I*sqrt(3)-3), 864/((I*sqrt(3)+3)^3*
   (I*sqrt(3)-3)^2), 1-I*sqrt(2), -(144*(sqrt(6)-(3*I)*sqrt(2)-I*sqrt(3)-
   3))/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1, 1, (-1+I*sqrt(3)+2^(2/3))
   /(I*sqrt(3)-3))+1/HeunG((2*I)*sqrt(3)/(I*sqrt(3)-3), 864/((I*sqrt(3)+3)^3*
   (I*sqrt(3)-3)^2), 1-I*sqrt(2), -(144*(sqrt(6)-(3*I)*sqrt(2)-I*sqrt(3)-
   3))/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1, 1, (-1+I*sqrt(3)+2^(2/3))
   /(I*sqrt(3)-3)))*x*HeunG((2*I)*sqrt(3)/(I*sqrt(3)-3), 864/((I*sqrt(3)+3)^3*
   (I*sqrt(3)-3)^2), 1-I*sqrt(2), -(144*(sqrt(6)-(3*I)*sqrt(2)-I*sqrt(3)-
   3))/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1, 1, (-1-2^(2/3)*x+I*sqrt(3)+2^(2
   /3))/(I*sqrt(3)-3))+_C2*x*HeunG((2*I)*sqrt(3)/(I*sqrt(3)-3), 
   864/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1+I*sqrt(2), (72*I)*(sqrt(6)-
   (3*I)*sqrt(2)-(2*I)*sqrt(3)-6)*sqrt(2)/((I*sqrt(3)+3)^3*(I*sqrt(3)-3)^2), 1, 
   1, (-1-2^(2/3)*x+I*sqrt(3)+2^(2/3))/(I*sqrt(3)-3))

By the way, why does Mathematice use these [FormalN] variables in output like this?

I do not know. I do not like them either. (most likely so they do not clash with your own symbols). MrWizard has a function here that can help remove them:

How to convert formal symbols to standard ones in an expression?

Update: After posting this, I asked Maple experts on the constant _C2 showing up in the solution of dsolve. Thanks to their help, it turned out the ode has wrong initial conditions. At $x=0$ the ode does not exist (the coefficients of the y' and y'' become zero). So the ODE generated by Mathematica, or the initial conditions that came with it, or both, are not correct. i.e. they lead to singularity at x=0.

This might indicate a bug in FindGeneratingFunction but I am not sure. Either way, the ODE given above can't be solved uniquely given the initial conditions that come with it.