Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
SAGE code for computing flow matrix for exp(z)-1
#3
(08/19/2009, 09:17 PM)jaydfox Wrote: Alas, it seems that Maxima and/or Lisp seems to crash for a large size, or after repeated calculations, due to small, hard-coded memory limit in the Maxima installation. According to discussions on the sage-support google group, it might require a recompilation of Maxima!

Well, I know ways to work around this, as there are other ways to compute the taylor series of composed functions, such as through an alternative library like gp/pari, but I don't have time to work on a workaround today. Sorry, folks!
Okay, so I decided to just calculate the taylor series of the iterations "by hand", using matrix math. It turns out to be way faster. I am now getting times comparable to what Andrew was getting in Mathematica. Hooray!

Code:
# Define the custom IDM:
def parabolic_idm(field, expr, x, size):
    mulm = matrix(field, size+1, size+1);
    abel = matrix(field, size+1, size+1);
    ret = matrix(field, size+1, size+1);
    cof = taylor(expr, x, 0, size).coeffs();
    for rrow in xrange(size+1):
        for ccol in xrange(len(cof)):
            col2 = cof[ccol][1]+rrow;
            if col2 < size+1:
                mulm[rrow, col2] = field(cof[ccol][0]);
    vec = vector(field, size+1);
    vec[0] = field(1);
    abel[0] = vec;
    for ii in xrange(size):
        vec = vec * mulm;
        abel[ii+1] = vec;
    vec = vector(field, size+1);
    vec[1] = field(1);
    ret[0] = vec;
    for ii in xrange(size):
        vec = vec * abel;
        ret[ii+1] = vec;
    return ret;

# Create a Lagrange Interpolation polynomial for each
# coefficient of the iteration function. Note that
# multiplying on the left by a row-vector of powers
# of n (where n is the desired iteration) will give
# a row vector which represents the coefficients of
# a power series in z. Alternatively, multiplying
# on the right by a column vector of powers of z
# will give the coefficients of a power series in n,
# which allows direct calculation of the nth iterate,
# centered at z.
def LagrangePolynomialMatrix(pidm, size):
    XI = matrix(ZZ, size+1);
    XI[0,0] = factorial(size);
    for ccol in xrange(1, size+1):
        XI[0, ccol] = (1-ccol)/ccol * XI[0, ccol-1]
    for rrow in xrange(1,size+1):
        for ccol in xrange(rrow, size+1):
            XI[rrow, ccol] =  (XI[rrow-1, ccol-1] + (1-ccol) * XI[rrow, ccol-1])/ccol
    PI = matrix(ZZ, size+1);
    PI[0,0] = 1;
    for rrow in xrange(1, size+1):
        for ccol in xrange(rrow+1):
            PI[rrow, ccol] = PI[rrow-1, ccol-1] - PI[rrow-1, ccol];
    return (XI*PI)*matrix(QQ, pidm)/factorial(size);

# test data: change or delete as needed
size = 100
time mtx = parabolic_idm(QQ, exp(x)-1, x, size)
time fma = LagrangePolynomialMatrix(mtx, size)

Update: It occurs to me that memory can be saved by eliminating the "mulm" variable and repurposing the "ret" variable. Thus:
Code:
# Define the custom IDM:
def parabolic_idm(field, expr, x, size):
    abel = matrix(field, size+1, size+1);
    ret = matrix(field, size+1, size+1);
    cof = taylor(expr, x, 0, size).coeffs();
    for rrow in xrange(size+1):
        for ccol in xrange(len(cof)):
            col2 = cof[ccol][1]+rrow;
            if col2 < size+1:
                ret[rrow, col2] = field(cof[ccol][0]);
    vec = vector(field, size+1);
    vec[0] = field(1);
    abel[0] = vec;
    for ii in xrange(size):
        vec = vec * ret;
        abel[ii+1] = vec;
    vec = vector(field, size+1);
    vec[1] = field(1);
    ret[0] = vec;
    for ii in xrange(size):
        vec = vec * abel;
        ret[ii+1] = vec;
    return ret;
~ Jay Daniel Fox
Reply


Messages In This Thread
RE: SAGE code for computing flow matrix for exp(z)-1 - by jaydfox - 08/19/2009, 11:47 PM

Possibly Related Threads...
Thread Author Replies Views Last Post
  Revisting my accelerated slog solution using Abel matrix inversion jaydfox 21 7,486 02/09/2019, 02:25 PM
Last Post: sheldonison
  C++ code for tet, ate and hexp MorgothV8 0 2,897 07/10/2014, 04:24 PM
Last Post: MorgothV8
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 6,515 08/15/2011, 03:14 PM
Last Post: Gottfried
  An incremental method to compute (Abel) matrix inverses bo198214 3 8,976 07/20/2010, 12:13 PM
Last Post: Gottfried
  Sage Question? rsgerard 1 3,894 05/09/2010, 11:40 AM
Last Post: bo198214
  Single-exp series computation code mike3 0 2,874 04/20/2010, 08:59 PM
Last Post: mike3
  computing teh last digits without computing the number deepinlife 3 5,622 02/24/2009, 09:09 AM
Last Post: deepinlife
  How to free memory in SAGE? jaydfox 2 5,372 12/21/2007, 06:38 PM
Last Post: andydude
  Convergence of matrix solution for base e jaydfox 6 8,636 12/18/2007, 12:14 AM
Last Post: jaydfox
  Computing Abel function at a given center jaydfox 10 11,945 11/30/2007, 06:44 PM
Last Post: andydude



Users browsing this thread: 1 Guest(s)