Ask Your Question

How to compute this exponential generating function?

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

Peter Luschny gravatar image

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);

This gives the requested triangle:

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

2 Answers

Sort by ยป oldest newest most voted

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

kcrisman gravatar image

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:

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 more


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

Peter Luschny gravatar imagePeter 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 gravatar imagerws ( 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 gravatar imagekcrisman ( 2014-10-31 21:02:46 +0200 )edit

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

Peter Luschny gravatar image

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.

edit flag offensive delete link 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.

kcrisman gravatar imagekcrisman ( 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

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

For what it's worth I shared an IP-worksheet at SageMath-Cloud:

Peter Luschny gravatar imagePeter 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 gravatar imageFrancis 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 gravatar imagePeter Luschny ( 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

Question Tools

1 follower


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

Seen: 1,196 times

Last updated: Oct 31 '14