# 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],, 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

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

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

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 - Yeah, I noticed that too. And it cycles through. I didn't have time to investigate further.

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)


[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

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.

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

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

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.