Ask Your Question

# How to compute this exponential generating function?

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 close merge delete

## 2 Answers

Sort by ยป oldest newest most voted

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.

more

## Comments

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

( 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()?

( 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.

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

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.

more

## 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.

( 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

( 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

( 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.

( 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.

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

## Your Answer

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

Add Answer

## Stats

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

Seen: 1,196 times

Last updated: Oct 31 '14