Finding recursive relations in Mathematica

[mnb96]

Hello,
I have a sequence of polynomials defined in the following way:
P_k(x) = \frac{\partial^k}{\partial x^k}e^{s(x)}\vert_{x=0}

Essentially the polynomial P_k is the k-th derivative of \exp(s(x)) evaluated at x=0. The function s(x) is a polynomial of 2nd degree in x.

In mathematica I define the polynomials with the following code:

D[Exp[s[x]], {x, k}]

For k=1..n one obtains a list of n polynomials P₁(x),...,P_n(x).
My question is: is it possible to ask Mathematica to find a recurrent relation that expresses any P_k as a function of the previous polynomials P_k-1, P_k-2... ? (assuming it exists).

A similar well-known problem exists for http://en.wikipedia.org/wiki/Hermite_polynomials#Definition"

Thanks.

 

[Simon_Tyler]

Try

In[1]:= s[x_] := a x^2+b x+c

In[2]:= poly[k_] = Simplify[SeriesCoefficient[Exp[s[x]], {x, 0, k}], k>0]

Out[2]= DifferenceRoot[Function[{y,n}, 
       {-2 a y[n] - b y[1+n] + (2+n) y[2+n]==0, y[0]==E^c,  y[1]==b E^c} ]]

The above shows the difference equation (recurrence relation) and the initial conditions that the polynomials must satisfy for any given quadratic s[x]. We can now use the recurrence relation to calculate the sequence. Below we compare calculating using poly[n] (the DifferenceRoot object) and Recurrence table

In[11]:= dr = Table[poly[n], {n, 0, 20}]; // Timing
         drExp = Expand[dr]; // Timing

Out[11]= {0.21, Null}
Out[12]= {1.57, Null}

In[13]:= rt = RecurrenceTable[{-2 a y[n] - b y[1 + n] + (2 + n) y[2 + n] == 0, 
               y[0] == E^c, y[1] == b E^c}, y, {n, 0, 20}]; // Timing
         rtExp = Expand[rt]; // Timing

Out[13]= {0.01, Null}
Out[14]= {1.66, Null}

In[15]:= dr == rt
Out[15]= True

In[16]:= drExp == rtExp
Out[16]= True

The problem is that both evaluating the DifferenceRoot poly[k] and using RecurrenceTable, whilst great for purely numeric functions, are horrible for symbolic functions such as the one above. the reason for this is that they don't simplify y[n] after every step and end up with a terrible mess, e.g.

In[17]:= dr[[6]]
         drExp[[6]]
Out[17]= 1/5 (2/3 a (2 a b E^c+1/2 b (2 a E^c+b^2 E^c))+1/4 b (a (2 a E^c+b^2 E^c)+1/3 b (2 a b E^c+1/2 b (2 a E^c+b^2 E^c))))
Out[18]= 1/2 a^2 b E^c+1/6 a b^3 E^c+(b^5 E^c)/120

Finally, compare calculating by taking derivatives with calculating using http://reference.wolfram.com/mathematica/tutorial/FunctionsThatRememberValuesTheyHaveFound.html" (in the latter, it's the Expand that makes the difference compared to the DifferenceRoot and RecurrenceTable options given above)

In[25]:= deriv=Table[1/k! D[Exp[s[x]], {x,k}] /. x->0, {k, 0, 100}]//Expand;//Timing
Out[25]= {0.36, Null}

In[26]:= Clear[y]
         y[n_Integer?Positive] := y[n] = 1/n (2 a y[-2+n] + b y[-1+n]) //Expand
         y[0] = E^c;
         y[1] = b E^c;

In[30]:= mem = Table[y[n], {n, 0, 100}]//Expand;//Timing
Out[30]= {0.3, Null}

In[31]:= deriv==mem
Out[31]= True

 

[Simon_Tyler]

Also, if you want the recurrence relation for the polynomials not evaluated at x=0, try

s[x_] := a x^2 + b x + c
poly[k_] = Simplify[SeriesCoefficient[Exp[s[x]], {x, x, k}], k > 0]

which returns

DifferenceRoot[Function[{y,n},
   {-2 a y[n] - (b + 2 a x) y[1+n] + (2+n) y[2+n]==0,
   y[0]==E^(c+b x+a x^2),  y[1]==(b+2 a x) E^(c+b x+a x^2)}
   ]][k]

 

[Simon_Tyler]

I must have had my brain turned off this morning, since I didn't mention that what you supplied was an exponential generating function. Just google "generating function recursion relation" to find stuff on how to approach these problems more generally. In particular the free book "generatingfunctionology" is worth reading (not that I can remember much of it)

 

[mnb96]

Thanks a lot for the answers!
Your hint about generating functions was also very useful. I hadn't thought about that!