How to compute this exponential generating function?

generating_function

asked 2014-10-30 19:50:23 +0200

Peter Luschny
548 ●13 ●35 ●56

Yesterday I struggled with the Narayana polynomials which can be computed with Sage (and the help of kcrisman) as

P = lambda n: hypergeometric([-n,-n+1],[2], 1/x).simplify_hypergeometric()
[expand(x^k*P(k)) for k in (0..7)]

Today I want to generate the coefficients of these polynomials by an exponential generating function. With Maple I have no problems.

egf := 1 + int((sqrt(t)*exp((1+t)*x)*BesselI(1,2*sqrt(t)*x))/x,x);
s := n -> n!*coeff(series(egf,x,n+2),x,n);
seq(print(seq(coeff(s(n),t,j),j=0..n)),n=0..6);

This gives the requested triangle:

1
0, 1
0, 1, 1
0, 1, 3, 1
0, 1, 6, 6, 1
0, 1, 10, 20, 10, 1
0, 1, 15, 50, 50, 15, 1

With Sage I tried:

from sage.symbolic.integration.integral import indefinite_integral
t = var('t')
h = lambda x, t: sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
egf = 1 + indefinite_integral(h(x,t),x)
taylor(egf, x, 0, 6)

This gives me a TypeError. Perhaps this is the wrong Ansatz.

Any help is appreciated.

edit retag flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2014-10-30 20:33:40 +0200

kcrisman
12222 ●41 ●135 ●254

updated 2014-10-31 14:34:27 +0200

Probably there is more to it than this, but why not just do

from sage.symbolic.integration.integral import indefinite_integral
h(x,t) = sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
egf = 1 + indefinite_integral(h,x)
taylor(egf, x, 0, 6)

This may not do it for you, but hopefully removes the error in question. Lambdas may not play well with integration.

Edit: probably the error was caused by using the indefinite integral directly.

sage: h(x,t) = sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
sage: egf = integrate(h,x) + 1

works fine. However,

sage: egf
(x, t) |--> sqrt(t)*integrate(bessel_I(1, 2*sqrt(t)*x)*e^((t + 1)*x), x) + 1

so Sage doesn't know how to do much with that, as the error message

sage: taylor(egf, x, 0, 6)
TypeError: ECL says: Error executing code in Maxima: taylor: unable to expand at a point specified in:
'integrate(bessel_i(1,2*sqrt(t)*x)*%e^((t+1)*x),x)

indicates. And the other way to get series (there is another, from Ginac) yields not much better:

sage: egf.series(x==0, 6)
AttributeError: 'MaximaLibElement' object has no attribute '_name'

since it doesn't know what to do with an unevaluated integral.

There may yet be a workaround, but I'm not conversant enough with this to say what it would be.

edit flag offensive delete link

Comments

TypeError: cannot coerce arguments: no canonical coercion from Callable function ring with arguments (x, t) to Symbolic Ring

Peter Luschny ( 2014-10-30 20:37:03 +0200 )edit

Funny: without the call to taylor(), I get instead AttributeError: 'sage.rings.integer.Integer' object has no attribute 'variables' from "egf.series(x==0,6)". Some global switch in taylor()?

rws ( 2014-10-31 17:00:53 +0200 )edit

@rws - Yeah, I noticed that too. And it cycles through. I didn't have time to investigate further.

kcrisman ( 2014-10-31 21:02:46 +0200 )edit

add a comment

answered 2014-10-31 19:20:25 +0200

Peter Luschny
548 ●13 ●35 ●56

This is not a solution but only a very awkward attempt which illustrates what is missing (coefficient_list is home-brewed).

egf = 1 + sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
tay = integral(taylor(egf, x, 0, 6),x)
r = coefficient_list(tay, x, true)
red = lambda i,o: o/i if i>0 else o
P = [red(ind, obj) for ind, obj in enumerate(r[1::])]  
for narayana_poly in P:
    print coefficient_list(narayana_poly, t)  

[1]
[0, 1]
[0, 1, 1]
[0, 1, 3, 1]
[0, 1, 6, 6, 1]
[0, 1, 10, 20, 10, 1]
[0, 1, 15, 50, 50, 15, 1]

I suspect Sage will not get many friends at the OEIS if it cannot cope with such elementary maneuvers with generating functions.

edit flag offensive delete link

Comments

Possibly. Do Taylor polys and integrals *always* commute - don't you need some sort of convergence criterion? If the homebrew is useful that would be good as well.

kcrisman ( 2014-10-31 21:09:39 +0200 )edit

Good question. On the other hand the generating functions are formal power series and the 'taylor' operation is in this context also to be understood in the formal sense. This is also the way I understand https://en.wikipedia.org/wiki/Generating_function

Peter Luschny ( 2014-11-01 11:07:55 +0200 )edit

For what it's worth I shared an IP-worksheet at SageMath-Cloud: https://cloud.sagemath.com/projects/b9c40d1b-4c81-4420-b74c-dc9b2c0318c6/files/OgfAndEgf.ipynb

Peter Luschny ( 2014-11-01 11:09:27 +0200 )edit

Probably worth noting that print oeis('Narayana')[0].programs()[0] gives a very simple formula for the coefficients.

Francis Clarke ( 2014-11-02 10:47:51 +0200 )edit

Francis, (1) the example was chosen to show that generating functions based on taylor cannot cope with integrals. (2) The 'very simple formula' Pari gives is unfortunately wrong if you use the definition given above in the OP.

Peter Luschny ( 2014-11-03 11:00:52 +0200 )edit

see more comments

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

How to compute this exponential generating function?

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

How to compute this exponential generating function? edit

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

How to compute this exponential generating function?