3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Rather than verifying just one instance of sf you can verify all of d1 with d1 == sf /@ Range[Length[d1]] $\endgroup$ – Bob Hanlon Nov 8 '14 at 5:22
4
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

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.

$\endgroup$
  • $\begingroup$ Thanks! I also used FriCAS' guessADE() to guess an algebraic differential equation for the generating function of my sequence and obtained {(x^3 - 3 x^2 + 3 x - 3) y''[x] + (3 x^2 - 6 x + 3) y'[x] + (3 x - 3) y[x] == 0, y[0] == 1, y'[0] == 1}, which Mma solves in terms of ordinary Gauss hypergeometric functions Hypergeometric2F1. $\endgroup$ – Eckhard Nov 8 '14 at 10:03
  • $\begingroup$ @Eckhard Ok. The ode you obtained from FriCAS looks ok. no singularity at x=0. This probably mean there might be a bug in Mathematica FindGeneratingFunction here. This might be a good example to send to support@wolfram.com to look at. $\endgroup$ – Nasser Nov 8 '14 at 10:11

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.