suggestions on how to do this

Discussion in 'Python' started by chris, Apr 27, 2005.

  1. chris

    chris Guest

    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
     
    chris, Apr 27, 2005
    #1
    1. Advertising

  2. chris

    Terry Reedy Guest

    "chris" <> wrote in message
    news:xNKbe.30981$...
    > 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
     
    Terry Reedy, Apr 27, 2005
    #2
    1. Advertising

  3. On Wed, 27 Apr 2005 11:34:53 GMT, "chris" <> wrote:

    >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
     
    Bengt Richter, Apr 27, 2005
    #3
  4. chris

    Kay Schluehr Guest

    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
     
    Kay Schluehr, Apr 27, 2005
    #4
  5. chris

    Guest

    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.
     
    , Apr 28, 2005
    #5
  6. chris

    Guest

    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
     
    , Apr 28, 2005
    #6
  7. chris

    bgs Guest

    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
     
    bgs, Apr 28, 2005
    #7
  8. chris

    bwaha Guest

    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
     
    bwaha, Apr 28, 2005
    #8
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. GriffithsJ

    Book suggestions

    GriffithsJ, Aug 15, 2003, in forum: ASP .Net
    Replies:
    1
    Views:
    903
    Wayne Wengert
    Aug 16, 2003
  2. Jeff Gaines

    Book Suggestions Please!

    Jeff Gaines, Sep 8, 2003, in forum: ASP .Net
    Replies:
    9
    Views:
    527
    William Ryan
    Sep 28, 2003
  3. JerryK
    Replies:
    2
    Views:
    600
    Ray Cassick \(Home\)
    Jan 29, 2004
  4. Victor Hannak
    Replies:
    1
    Views:
    567
    Mike Treseler
    Nov 25, 2003
  5. Maxim
    Replies:
    0
    Views:
    420
    Maxim
    Jul 7, 2003
Loading...

Share This Page