suggestions on how to do this

[chris]

The problem I have is as follows:

I have a recursive function b(k)

b(k) = -(A/k**2)*(b(k-2) - b(k-5))
k<0, b(k)=0
k=0, b(k)=1
k=1, b(k)=0

eg. b(2) = -A/4
b(3) = 0
b(4) = A**2/64

note that as k increases b(k) can itself be a sum of terms in powers of A
rather than a single power of A in the examples above.

Summing all terms and equating to zero gives:

F= sum b(k) = 0 for all k = 0, infinity

When this is expanded I get a polynomial F(A). I want to determine the
coefficients of the polynomial so that I can find the roots of the function
F up to a specified order of A.


I have yet to code this but I was hoping for some ideas on how to do this
reasonably.

I figure I can compute each b(k) and store the numeric value(s) and
associated powers of A. Then collect coefficients for like powers of A.
Finally I have a set of polynomial coefficients in A which I can pass to
scipy.base.roots()

Any suggestions on how I might do this efficiently? I have no doubt I can
get this done with brute force, but I would prefer to explore more elegant
means which I look to the masters for.

tia

 

[Terry Reedy]

I have a recursive function b(k)

b(k) = -(A/k**2)*(b(k-2) - b(k-5))

This is specifically called a recurrence relation/

[snip]

When this is expanded I get a polynomial F(A). I want to determine the
coefficients of the polynomial so that I can find the roots of the
function
F up to a specified order of A.

This is a math/numerical analysis problem. If you don't get an answer
here, I'd look for a text of recurrence relations and polynomials, and look
for a math/numerical analysis newsgroup. sci.mathematics?

Note: your subject line is so vague that a specialist who can tell you more
could but only selectively opens threads could easily miss this. It is
only happenstance that I didn't.

Terry J. Reedy

 

[Bengt Richter]

The problem I have is as follows:

I have a recursive function b(k)

b(k) = -(A/k**2)*(b(k-2) - b(k-5))
k<0, b(k)=0
k=0, b(k)=1
k=1, b(k)=0

eg. b(2) = -A/4
b(3) = 0
b(4) = A**2/64

note that as k increases b(k) can itself be a sum of terms in powers of A
rather than a single power of A in the examples above.

Summing all terms and equating to zero gives:

F= sum b(k) = 0 for all k = 0, infinity

When this is expanded I get a polynomial F(A). I want to determine the
coefficients of the polynomial so that I can find the roots of the function
F up to a specified order of A.


I have yet to code this but I was hoping for some ideas on how to do this
reasonably.

I figure I can compute each b(k) and store the numeric value(s) and
associated powers of A. Then collect coefficients for like powers of A.
Finally I have a set of polynomial coefficients in A which I can pass to
scipy.base.roots()

Any suggestions on how I might do this efficiently? I have no doubt I can
get this done with brute force, but I would prefer to explore more elegant
means which I look to the masters for.

Does this look right?

b(-5) -> 0
b(-4) -> 0
b(-3) -> 0
b(-2) -> 0
b(-1) -> 0
b(0) -> 0
b(1) -> 0
b(2) -> -A/4
b(3) -> 0
b(4) -> A**2/64
b(5) -> A/25
b(6) -> -A**3/2304
b(7) -> -29*A**2/4900
b(8) -> A**4/147456
b(9) -> 563*A**3/2116800
b(10) -> A**2/2500 -A**5/14745600
b(11) -> -5927*A**4/1024531200
b(12) -> -43*A**3/980000 +A**6/2123366400
b(13) -> 824003*A**5/11081329459200
b(14) -> 16397*A**4/10372320000 -A**7/416179814400
b(15) -> A**3/562500 -1260403*A**6/1994639302656000

Regards,
Bengt Richter

 

[Kay Schluehr]

Seems You are looking for a generating function of a recurrence
relation. There is a standard text about this topic written by Herbert
S. Wilf downloadable from his homapage:

http://www.cis.upenn.edu/~wilf/

Regards,
Kay

 

[cfriedl]

Hi Terry

I apprecaite the advice. Briefly I'm familiar with the math (its an
eigenvalue problem in fluid flow) but not so much with python (3 months
on and off), hence my post here. I'm looking for python advice on how
to implement this effectively. I love python and would like to use it
well.

As for vague subject, I've not found a strong correlation between
subject line and number or depth of responses. Perhaps today the
curious looked at it. But at other times, even with a quite specific
subject, no replies come. Such is life.

 

[cfriedl]

Hi Bengt

OK, that's right. So I'm curious how you did this. And any comments on
how to collect coefficients of like powers in your result.

Be Well and Happy Always

Chris

 

[bgs]

You could probably use scipy.base.polynomial, but it's easy enough to
implement a polynomial yourself. Just use a dict-- each key represents
the power and each value the coefficient of the polynomial.

You didn't say exactly how efficient you need this. It takes only a
couple seconds to sum 100 of the b(k)'s using the implementation below.
This gets you roots out to about A=500.

Perhaps there are bugs in the below implementation. Either way, you
can compare your results and its results and if they're not the same,
well then there's a problem.


------------------------------------------
from rational import rational #get from bitconjurer.org/rational.py

def add(a, b, s=1):
c = {}
for k, v in b.items():
if not a.has_key(k):
c[k] = s*v
for k, v in a.items():
vb = b.get(k, 0)
c[k] = a[k] + s*vb
return c

def raise1(a):
b = {}
for k, v in a.items():
b[k+1] = v
return b

def scale(a, s):
b = {}
for k, v in a.items():
b[k] = v*s
return b

def eval(x, p):
s = 0.0
for k, v in p.items():
s = s + float(v.num)/float(v.den)*x**k
return s

b = {-3 : {}, -2 : {}, -1 : {}, 0 : {0:rational(1,1)}, 1 : {}}

N = 100
for k in range(2,N):
b[k] = scale(raise1(add(b[k-5],b[k-2],-1)),rational(1,k**2))
o = [b[0]]
for k in range(1, N):
o.append(add(o[-1], b[k]))

for x in range(0,800):
print x, eval(x, o[-3]), eval(x, o[-2]), eval(x, o[-1])
# o[-1],o[-2], and o[-3] start to split around x = 500

 

[bwaha]

Hi bgs

Many thanks. This generates the correct coefficients. I am studying
your implementation now. I've not used a dictionary of dictionaries
before so there's a bit of a learning curve going on right now. However
I can see that b[k] holds the relevant info (coefficients and powers)
so I can easily modify it to collect terms with like powers. As for
speed, its not a concern. Thanks again.

bwaha Chris