ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 03 Nov 2014 08:09:19 -0600How to compute this exponential generating function?http://ask.sagemath.org/question/24713/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.Thu, 30 Oct 2014 13:50:23 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/Answer by Peter Luschny for <p>Yesterday I struggled with the Narayana polynomials which can be computed with Sage (and the help of kcrisman) as</p>
<pre><code>P = lambda n: hypergeometric([-n,-n+1],[2], 1/x).simplify_hypergeometric()
[expand(x^k*P(k)) for k in (0..7)]
</code></pre>
<p>Today I want to generate the coefficients of these polynomials by an exponential generating function.
With Maple I have no problems.</p>
<pre><code>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);
</code></pre>
<p>This gives the requested triangle:</p>
<pre><code>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
</code></pre>
<p>With Sage I tried:</p>
<pre><code>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)
</code></pre>
<p>This gives me a TypeError. Perhaps this is the wrong Ansatz. </p>
<p>Any help is appreciated.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?answer=24730#post-id-24730This 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](https://oeis.org/A090181) if it cannot cope with such elementary maneuvers with generating functions.Fri, 31 Oct 2014 13:20:25 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?answer=24730#post-id-24730Comment by Peter Luschny for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24737#post-id-24737Good 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_functionSat, 01 Nov 2014 05:07:55 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24737#post-id-24737Comment by kcrisman for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24758#post-id-24758PL - yes, I would agree, except that in this case taylor isn't generating formal power series necessarily. In particular, I don't think we have a generating function class different from power series. In any case, that is probably not really germane to the discussion, which you correctly identify as taylor series and integrals not playing well together in general, perhaps especially when there are things like Bessel functions involved that are (presumably) nontrivial to symbolically integrate.Mon, 03 Nov 2014 08:09:19 -0600http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24758#post-id-24758Comment by Peter Luschny for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24748#post-id-24748Francis, (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.Mon, 03 Nov 2014 04:00:52 -0600http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24748#post-id-24748Comment by kcrisman for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24734#post-id-24734Possibly. 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.Fri, 31 Oct 2014 15:09:39 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24734#post-id-24734Comment by Peter Luschny for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24738#post-id-24738For what it's worth I shared an IP-worksheet at SageMath-Cloud: https://cloud.sagemath.com/projects/b9c40d1b-4c81-4420-b74c-dc9b2c0318c6/files/OgfAndEgf.ipynbSat, 01 Nov 2014 05:09:27 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24738#post-id-24738Comment by Francis Clarke for <p>This is not a solution but only a very awkward attempt which
illustrates what is missing (coefficient_list is home-brewed). </p>
<pre><code>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]
</code></pre>
<p>I suspect Sage will not get many friends at the <a href="https://oeis.org/A090181">OEIS</a> if it cannot cope with such elementary maneuvers with generating functions.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24742#post-id-24742Probably worth noting that print oeis('Narayana')[0].programs()[0] gives a very simple formula for the coefficients.Sun, 02 Nov 2014 03:47:51 -0600http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24742#post-id-24742Answer by kcrisman for <p>Yesterday I struggled with the Narayana polynomials which can be computed with Sage (and the help of kcrisman) as</p>
<pre><code>P = lambda n: hypergeometric([-n,-n+1],[2], 1/x).simplify_hypergeometric()
[expand(x^k*P(k)) for k in (0..7)]
</code></pre>
<p>Today I want to generate the coefficients of these polynomials by an exponential generating function.
With Maple I have no problems.</p>
<pre><code>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);
</code></pre>
<p>This gives the requested triangle:</p>
<pre><code>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
</code></pre>
<p>With Sage I tried:</p>
<pre><code>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)
</code></pre>
<p>This gives me a TypeError. Perhaps this is the wrong Ansatz. </p>
<p>Any help is appreciated.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?answer=24714#post-id-24714Probably 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.Thu, 30 Oct 2014 14:33:40 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?answer=24714#post-id-24714Comment by kcrisman for <p>Probably there is more to it than this, but why not just do </p>
<pre><code>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)
</code></pre>
<p>This may not do it for you, but hopefully removes the error in question. Lambdas may not play well with integration.</p>
<p>Edit: probably the error was caused by using the indefinite integral directly.</p>
<pre><code>sage: h(x,t) = sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
sage: egf = integrate(h,x) + 1
</code></pre>
<p>works fine. However, </p>
<pre><code>sage: egf
(x, t) |--> sqrt(t)*integrate(bessel_I(1, 2*sqrt(t)*x)*e^((t + 1)*x), x) + 1
</code></pre>
<p>so Sage doesn't know how to do much with that, as the error message </p>
<pre><code>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)
</code></pre>
<p>indicates. And the other way to get series (there is another, from Ginac) yields not much better:</p>
<pre><code>sage: egf.series(x==0, 6)
AttributeError: 'MaximaLibElement' object has no attribute '_name'
</code></pre>
<p>since it doesn't know what to do with an unevaluated integral.</p>
<p>There may yet be a workaround, but I'm not conversant enough with this to say what it would be.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24733#post-id-24733@rws - Yeah, I noticed that too. And it cycles through. I didn't have time to investigate further.Fri, 31 Oct 2014 15:02:46 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24733#post-id-24733Comment by rws for <p>Probably there is more to it than this, but why not just do </p>
<pre><code>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)
</code></pre>
<p>This may not do it for you, but hopefully removes the error in question. Lambdas may not play well with integration.</p>
<p>Edit: probably the error was caused by using the indefinite integral directly.</p>
<pre><code>sage: h(x,t) = sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
sage: egf = integrate(h,x) + 1
</code></pre>
<p>works fine. However, </p>
<pre><code>sage: egf
(x, t) |--> sqrt(t)*integrate(bessel_I(1, 2*sqrt(t)*x)*e^((t + 1)*x), x) + 1
</code></pre>
<p>so Sage doesn't know how to do much with that, as the error message </p>
<pre><code>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)
</code></pre>
<p>indicates. And the other way to get series (there is another, from Ginac) yields not much better:</p>
<pre><code>sage: egf.series(x==0, 6)
AttributeError: 'MaximaLibElement' object has no attribute '_name'
</code></pre>
<p>since it doesn't know what to do with an unevaluated integral.</p>
<p>There may yet be a workaround, but I'm not conversant enough with this to say what it would be.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24728#post-id-24728Funny: 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()?Fri, 31 Oct 2014 11:00:53 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24728#post-id-24728Comment by Peter Luschny for <p>Probably there is more to it than this, but why not just do </p>
<pre><code>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)
</code></pre>
<p>This may not do it for you, but hopefully removes the error in question. Lambdas may not play well with integration.</p>
<p>Edit: probably the error was caused by using the indefinite integral directly.</p>
<pre><code>sage: h(x,t) = sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
sage: egf = integrate(h,x) + 1
</code></pre>
<p>works fine. However, </p>
<pre><code>sage: egf
(x, t) |--> sqrt(t)*integrate(bessel_I(1, 2*sqrt(t)*x)*e^((t + 1)*x), x) + 1
</code></pre>
<p>so Sage doesn't know how to do much with that, as the error message </p>
<pre><code>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)
</code></pre>
<p>indicates. And the other way to get series (there is another, from Ginac) yields not much better:</p>
<pre><code>sage: egf.series(x==0, 6)
AttributeError: 'MaximaLibElement' object has no attribute '_name'
</code></pre>
<p>since it doesn't know what to do with an unevaluated integral.</p>
<p>There may yet be a workaround, but I'm not conversant enough with this to say what it would be.</p>
http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24715#post-id-24715TypeError: cannot coerce arguments: no canonical coercion from Callable function ring with arguments (x, t) to Symbolic RingThu, 30 Oct 2014 14:37:03 -0500http://ask.sagemath.org/question/24713/how-to-compute-this-exponential-generating-function/?comment=24715#post-id-24715