Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Finding recursive relations in Mathematica

  1. Aug 13, 2011 #1
    Hello,
    I have a sequence of polynomials defined in the following way:
    [tex]P_k(x) = \frac{\partial^k}{\partial x^k}e^{s(x)}\vert_{x=0}[/tex]

    Essentially the polynomial Pk is the k-th derivative of [itex]\exp(s(x))[/itex] 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 P1(x),...,Pn(x).
    My question is: is it possible to ask Mathematica to find a recurrent relation that expresses any Pk as a function of the previous polynomials Pk-1, Pk-2... ? (assuming it exists).

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

    Thanks.
     
    Last edited by a moderator: Apr 26, 2017
  2. jcsd
  3. Aug 13, 2011 #2
    Try

    Code (Text):
    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

    Code (Text):

    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.

    Code (Text):

    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" [Broken] (in the latter, it's the Expand that makes the difference compared to the DifferenceRoot and RecurrenceTable options given above)

    Code (Text):

    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
     
    Last edited by a moderator: May 5, 2017
  4. Aug 13, 2011 #3
    Also, if you want the recurrence relation for the polynomials not evaluated at x=0, try

    Code (Text):

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

    Code (Text):
    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]
     
  5. Aug 14, 2011 #4
    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)
     
  6. Aug 16, 2011 #5
    Thanks a lot for the answers!
    Your hint about generating functions was also very useful. I hadn't thought about that!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Finding recursive relations in Mathematica
Loading...