Laplace Inverse Transform Implementation

Laplace Inverse Transform Implementation - Printable Version

+- Python Forum (https://python-forum.io)
+-- Forum: Python Coding (https://python-forum.io/forum-7.html)
+--- Forum: General Coding Help (https://python-forum.io/forum-8.html)
+--- Thread: Laplace Inverse Transform Implementation (/thread-927.html)

Pages: 1 2

Laplace Inverse Transform Implementation - jafar1363 - Nov-15-2016

Hello

My name is Jafar. In a part of my research I need to use DE HOOG inverse Laplace transform algorithm in Python.

I found a code from this link:

Code

#Copyright (C) 2009  Karl Bandilla
#
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program.  If not, see http://www.gnu.org/licenses/.

from numpy import *

def HoogTransform(t, gamma, bigT, N, FUN, meth, init, d, work):
    '''Numerrical Laplace inversion using the Hoog method.
       t     (float variable) independent variable at which inverse
             function is to be evaluated
       gamma (float variable) parameter of inversion integral, governing
             accuracy of inversion
       bigT  (float variable) parameter used to discretize inversion
             integral, thus governing discretization error
       N     (integer variable) must be >= 2 and should be even; number
             of terms to be used in sum for approximation, thus
             governing truncation error
       FUN   name of user-supplied function which evaluates the complex
             values of transform fbar(p) for complex values of Laplace
             transform variable p
       meth  (integer variable) determines inversion method: 1=epsilon
             method; 2=ordinary quotient difference algorithm;
             3=acceleratedquotient difference algorithm
       init  (integer variable) =1 indicates this is the first call
             of the algorithm with these values of gamma, bigT, N and
             FUN. !=1 means array is not recomputed
       d     (complex array dimensioned in calling function) contains
             continued fraction coefficient if init >1 and meth=2 or 3
       work  (complex array dimensioned in calling function) work space
             used for calculating successive diagonals of the table

       This is a translation of the FORTRAN subroutine LAPADC originally by
       J.H. Knight. See code for more information.
       Requires numpy.'''
    #This is a translation of the FORTRAN subroutine LAPADC by J.H. Knight
    #of the Division of Environmental Mechanics, CSIRO, Canberra, ACT2601,
    #Australia originally coded in August 1986. The function performs a
    #numerical LaPlace inversion using the Hoog method. Two papers are
    #referenced in the original FORTRAN code:
    #DeHoog, F.R., Knight, J.H. and Stokes, A.N. An improved method for
    #numerical inversion of LaPlace transforms. SIAM J. SCI. STA. COMPUT.,
    # 3, PP. 357-366, 1982.
    #DeHoog, F.R., Knight, J.H. and Stokes, A.N. Subroutine LAPACC for LaPlace
    #inversion by acceleration of convergence of complex sum. ACM TRANS.
    #MATH. SOFTWARE, (manuscript in preperation).
    #Translated by Karl W. Bandilla in July 2009.
    #Comments are taken from the FORTRAN code

    #Following is the header of the original FORTRAN code

#      SUBROUTINE LAPADC(T,GAMMA,BIGT,N,F,FUN,METH,INIT,D,WORK)
#
#-----------------------------------------------------------------------
#
#     AUTHOR
#
#         J.H. KNIGHT
#         DIVISION OF ENVIRONMENTAL MECHANICS
#         CSIRO, CANBERRA, ACT 2601, AUSTRALIA
#         AUGUST, 1986
#
#     REFERENCES
#
#         DE HOOG, F.R., KNIGHT, J.H. AND STOKES, A.N.
#         AN IMPROVED METHOD FOR NUMERICAL INVERSION OF LAPLACE
#         TRANSFORMS. SIAM J. SCI. STAT. COMPUT., 3, PP. 357-366, 1982.
#
#         DE HOOG, F.R., KNIGHT, J.H. AND STOKES, A.N.
#         SUBROUTINE LAPACC FOR LAPLACE TRANSFORM INVERSION BY
#         ACCELERATION OF CONVERGENCE OF COMPLEX SUM.
#         ACM TRANS. MATH. SOFTWARE, (MANUSCRIPT IN PREPARATION).
#
#     ABSTRACT
#
#         LAPADC COMPUTES AN APPROXIMATION TO THE INVERSE LAPLACE
#         TRANSFORM F(T) CORRESPONDING TO THE LAPLACE TRANSFORM
#         FBAR(P), USING A RATIONAL APPROXIMATION TO THE FOURIER
#         SERIES RESULTING WHEN THE INVERSION INTEGRAL IS
#         DISCRETIZED USING THE TRAPEZOIDAL RULE.
#         THE RATIONAL APPROXIMATION IS A DIAGONAL PADE APPROXIMATION
#         COMPUTED USING EITHER THE EPSILON ALGORITHM, IN WHICH THE
#         EPSILON TABLE MUST BE RECOMPUTED AT EACH VALUE OF T, OR
#         USING THE QUOTIENT-DIFFERENCE ALGORITHM, IN WHICH THE
#         QUOTIENT-DIFFERENCE ALGORITHM IS USED TO COMPUTE THE
#         COEFFICIENTS OF THE ASSOCIATED CONTINUED FRACTION, AND THE
#         CONTINUED FRACTION IS EVALUATED AT EACH VALUE OF T.
#         IN ADDITION, THE EVALUATION OF THE CONTINUED FRACTION IS
#         MADE MORE ACCURATE WITH AN IMPROVED ESTIMATE OF THE
#         REMAINDER.
#
#     METHOD
#
#         A REAL VALUED INVERSE FUNCTION F(T) IS THE REAL PART OF THE
#         INTEGRAL
#                                INFINITY
#
#         (ONE/PI)*EXP(GAMMA*T)* INTEGRAL FBAR(GAMMA+I*W)*EXP(I*W*T)*DW
#
#                                W=ZERO
#
#         WITH FBAR THE LAPLACE TRANSFORM, GAMMA A REAL PARAMETER, AND
#         I**2=-1.  W IS THE VARIABLE OF INTEGRATION.
#         IF THIS IS DISCRETIZED USING THE TRAPEZOIDAL RULE WITH STEP
#         SIZE PI/BIGT, A FOURIER SERIES RESULTS, AND F(T) IS THE SUM
#         OF THE DISCRETIZATION ERROR
#
#               INFINITY
#
#               SUM     EXP(-TWO*GAMMA*K*BIGT)*F(T+TWO*K*BIGT)
#
#               K=1
#
#         AND THE REAL PART OF THE FOURIER SERIES
#
#                                    INFINITY
#
#           (ONE/BIGT)*EXP(GAMMA*T)* SUM  A(K)*EXP(I*K*PI*T/BIGT)
#
#                                    K=0
#
#         WITH THE COEFFICIENTS (A(K),K=0,INFINITY) GIVEN BY
#
#              A(0) = HALF*FBAR(GAMMA),
#
#              A(K) = FBAR(GAMMA+I*K*PI/BIGT),   K=1,INFINITY.
#
#         THIS EXPRESSION CAN ALSO BE CONSIDERED AS A POWER SERIES IN
#         THE VARIABLE Z=EXP(I*PI*T/BIGT), AND WE INTRODUCE A
#         TRUNCATION ERROR WHEN WE USE A FINITE POWER SERIES
#
#                        M2
#
#              V(Z,M2) = SUM  A(K)*Z**K,       M2 = M*2.
#
#                        K=0
#
#         THE EPSILON ALGORITHM CALCULATES THE DIAGONAL PADE
#         APPROXIMATION TO THIS FINITE POWER SERIES
#
#                        M                M
#
#              V(Z,M2) = SUM  B(K)*Z**K / SUM  C(K)*Z**K
#
#                        K=0              K=0
#
#                                    + HIGHER ORDER TERMS,   C(0)=ONE.
#
#         THE CALCULATION MUST BE DONE FOR EACH NEW VALUE OF T.
#         THIS METHOD IS USED FOR METH=1.
#
#         THE QUOTIENT-DIFFERENCE (QD) ALGORITHM IS ANOTHER WAY
#         OF CALCULATING THIS DIAGONAL PADE APPROXIMATION, IN THE FORM
#         OF A CONTINUED FRACTION
#
#         V(Z,M2) = D(0)/(ONE+D(1)*Z/(ONE+...+D(M2-1)*Z/(ONE+R2M(Z)))),
#
#         WHERE R2M(Z) IS THE UNKNOWN REMAINDER.
#         THE QD ALGORITHM GIVES THE COEFFICIENTS (D(K),K=0,M2), AND
#         THE CONTINUED FRACTION IS THEN EVALUATED BY FORWARD RECURRENCE,
#         FOR EACH NEW VALUE OF T, WITH R2M(Z) TAKEN AS D(M2)*Z.
#         THE ANSWERS SHOULD BE IDENTICAL TO THOSE GIVEN BY THE
#         EPSILON ALGORITHM, APART FROM ROUNDOFF ERROR.
#         THIS METHOD IS USED FOR METH=2.
#
#         IN ADDITION, THE CONVERGENCE OF THE CONTINUED FRACTION CAN
#         BE ACCELERATED WITH AN IMPROVED ESTIMATE OF THE REMAINDER,
#         SATISFYING
#
#               R2M + ONE + (D(M2-1)-D(M2))*Z = D(M2)*Z/R2M
#
#         THIS METHOD IS USED FOR METH=3.
#
#     USAGE
#
#         THE USER MUST CHOOSE THE PARAMETERS BIGT AND GAMMA WHICH
#        DETERMINE THE DISCRETIZATION ERROR, AND N(=M2) WHICH
#         DETERMINES THE TRUNCATION ERROR AND ROUNDOFF ERROR, GIVEN
#         THE CHOICE OF BIGT AND GAMMA.  AS N IS INCREASED, THE
#         TRUNCATION ERROR DECREASES BUT THE ROUNDOFF ERROR INCREASES,
#         SINCE THE EPSILON AND QD ALGORITHMS ARE NUMERICALLY UNSTABLE.
#         IT IS DESIRABLE TO USE DOUBLE PRECISION FOR THE CALCULATIONS,
#         UNLESS THE COMPUTER KEEPS AT LEAST 12 SIGNIFICANT DECIMAL
#         DIGITS IN SINGLE PRECISION.
#         METH=1 USES THE EPSILON ALGORITHM, METH=2 USES THE QUOTIENT-
#         DIFFERENCE ALGORITHM, AND METH=3 USES THE QD ALGORITHM WITH
#         IMPROVED ESTIMATE OF THE REMAINDER ON EVALUATION.
#
#         AN ESTIMATE OF THE OPTIMAL VALUE OF N FOR GIVEN BIGT AND
#         GAMMA CAN BE GOT BY VARYING N AND CHOOSING A VALUE SLIGHTLY
#         LESS THAN THAT WHICH CAUSES THE RESULTS FOR METH=1 TO VARY
#         SIGNIFICANTLY FROM THOSE FOR METH=2.
#
#         CALCULATIONS USING THE QD ALGORITHM ARE MOST EFFICIENT
#         AND MOST ACCURATE IF LAPADC IS CALLED WITH METH=3, AND WITH
#         INIT=1 ON THE FIRST CALL WITH NEW VALUES OF BIGT, GAMMA, N
#         OR FUN, AND WITH INIT=2 ON SUBSEQUENT CALLS WITH NEW VALUES
#         OF T AND THE OTHER PARAMETERS UNCHANGED.
#
#
#     DESCRIPTION OF ARGUMENTS
#
#         T      -  REAL (DOUBLE PRECISION) VARIABLE.  ON INPUT, VALUE
#                   OF INDEPENDENT VARIABLE AT WHICH INVERSE FUNCTION
#                   F(T) IS TO BE APPROXIMATED.  UNCHANGED ON OUTPUT.
#         GAMMA  -  REAL (DOUBLE PRECISION) VARIABLE.  ON INPUT,
#                   CONTAINS PARAMETER OF INVERSION INTEGRAL, GOVERNING
#                   ACCURACY OF INVERSION.  UNCHANGED ON OUTPUT.
#         BIGT   -  REAL (DOUBLE PRECISION) VARIABLE.  ON INPUT,
#                   CONTAINS PARAMETER USED TO DISCRETIZE INVERSION
#                   INTEGRAL. GOVERNS DISCRETIZATION ERROR.
#                   UNCHANGED ON OUTPUT.
#         N      -  INTEGER VARIABLE, SHOULD BE EVEN .GE. 2.  ON INPUT,
#                   CONTAINS NUMBER OF TERMS TO BE USED IN SUM FOR
#                   APPROXIMATION.  DETERMINES TRUNCATION ERROR.
#                   UNCHANGED ON OUTPUT.
#         F      -  REAL (DOUBLE PRECISION) VARIABLE.  ON OUTPUT,
#                   CONTAINS VALUE OF INVERSE FUNCTION F(T), EVALUATED
#                   AT T.
#         FUN    -  NAME OF USER-SUPPLIED SUBROUTINE WHICH EVALUATES
#                   (DOUBLE) COMPLEX VALUES OF TRANSFORM FBAR(P)
#                   FOR (DOUBLE) COMPLEX VALUES OF LAPLACE
#                   TRANSFORM VARIABLE P.  CALLED IF METH .EQ. 1 .OR.
#                   INIT .EQ. 1.  MUST BE DECLARED EXTERNAL IN
#                   (SUB)PROGRAM CALLING LAPADC.
#         METH   -  INTEGER VARIABLE, DETERMINING METHOD TO BE USED TO
#                   ACCELERATE CONVERGENCE OF SUM.
#                   ON INPUT, METH .EQ. 1 INDICATES THAT THE EPSILON
#                   ALGORITHM IS TO BE USED.
#                   METH .EQ. 2 INDICATES THAT THE ORDINARY QUOTIENT
#                   DIFFERENCE ALGORITHM IS TO BE USED.
#                   METH .EQ. 3 INDICATES THAT THE QUOTIENT DIFFERENCE
#                   ALGORITHM IS TO BE USED, WITH FURTHER ACCELERATION
#                   OF THE CONVERGENCE OF THE CONTINUED FRACTION.
#                   UNCHANGED ON OUTPUT.
#         INIT   -  INTEGER VARIABLE. ON INPUT, USED ONLY IF Q-D
#                   ALGORITHM IS SPECIFIED (METH .EQ. 2 OR 3).
#                   INIT .EQ. 1 INDICATES THAT THIS IS THE FIRST CALL
#                   OF THE ALGORITHM WITH THESE VALUES OF GAMMA, BIGT,
#                   N, FUN.  IN THIS CASE, THE VALUES OF THE CONTINUED
#                   FRACTION COEFFICIENTS D(0:N) ARE RECALCULATED, AND
#                   THE CONTINUED FRACTION IS EVALUATED.
#                   INIT .NE. 1 INDICATES THAT THE VALUES OF GAMMA,
#                   BIGT, N, AND FUN ARE UNCHANGED SINCE THE LAST CALL.
#                   IN THIS CASE, THE INPUT VALUES OF THE CONTINUED
#                   FRACTION COEFFICIENTS D(0:N) ARE USED, AND THE
#                   CONTINUED FRACTION IS EVALUATED.
#                   UNCHANGED ON OUTPUT.
#         D      -  (DOUBLE) COMPLEX ARRAY, OF VARIABLE DIMENSION (0:N).
#                   MUST BE DIMENSIONED IN CALLING PROGRAM.
#                   ON INPUT, IF METH .EQ. 2 OR 3 .AND. INIT .NE. 1,
#                   MUST CONTAIN CONTINUED FRACTION COEFFICIENTS
#                   OUTPUT ON A PREVIOUS CALL.
#                   ON OUTPUT, IF METH .EQ. 2 OR 3, CONTAINS CONTINUED
#                   FRACTION COEFFICIENTS CALCULATED ON THIS CALL
#                   (INIT .EQ. 1) OR AS SUPPLIED ON INPUT (INIT .NE. 1).
#         WORK   -  (DOUBLE) COMPLEX ARRAY OF VARIABLE DIMENSION (0:N).
#                   MUST BE DIMENSIONED IN CALLING PROGRAM.
#                   WORK SPACE USED FOR CALCULATING SUCCESSIVE DIAGONALS
#                   OF THE TABLE IF METH .EQ. 1 .OR. INIT .EQ. 1.
#                   ON OUTPUT, CONTAINS LAST CALCULATED DIAGONAL OF
#                   TABLE.
#
#     SUBROUTINE CALLED
#
#           FUN(P,FBAR) - USER-SUPPLIED SUBROUTINE WHICH EVALUATES
#                   (DOUBLE) COMPLEX VALUES OF TRANSFORM FBAR(P)
#                   FOR (DOUBLE) COMPLEX VALUES OF LAPLACE
#                   TRANSFORM VARIABLE P.  CALLED IF METH .EQ. 1 .OR.
#                   INIT .EQ. 1.  MUST BE DECLARED EXTERNAL IN
#                   (SUB)PROGRAM WHICH CALLS SUBROUTINE LAPADC.


    if N < 2:
        print 'Order too low (N at least 2)'
        return -1.0

    zeror = 0.0
    zero = 0.0 + 0.0j
    half = 0.5 + 0.0j
    one = 1.0 + 0.0j
#    work = zeros(N+1,dtype=complex128)
#    d = zeros(N+1,dtype=complex128)
    #make order multiple of 2
    M2 = (N / 2) * 2
    factor = pi / bigT
 #   print 'factor', factor
    argR = t * factor
    Z = cos(argR) + sin(argR) * 1j
    if meth == 1 or init == 1:
        #new: get all A's first
        ##########################
        #need to think of way to do this
        ###########################
        #calculate table if epsilon algorithm is to be used, or
        #if this is the first call to quotient-difference algorithm
        #with these parameters and inverse transform
        A = FUN(gamma + zeror * 1j)
        AOld = half * A
        if AOld == 0.0:
            return 0.0
        argI = factor
        A = FUN(gamma + argI * 1j)
        #initialize table entries
        if meth == 1:
            ztoj = Z
            term = A * ztoj
            work[0] = AOld + term
            work[1] = one / term
        else:
            d[0] = AOld
            work[0] = zero
            work[1] = A / AOld
            d[1] = -1.0 * work[1]
            AOld = A
        #calculate successive diagonals of table
        for j in range(2, N + 1):
            #initialize calculation of diagonal
            old2 = work[0]
            old1 = work[1]
            argI += factor
            A = FUN(gamma + argI * 1j)
         #   print 'A', A, A - AOld
            if A == 0:
                return 0.0
            if meth == 1:
                #calculate next term and sum of power series
                ztoj = Z * ztoj
                term = A * ztoj
                work[0] += term
                work[1] = one / term
            else:
                work[0] = zero
                work[1] = A / AOld
                AOld = A
        #    print work[0:10]
            #calculate diagonal using Rhombus rules
            for i in range(2, j+1):
                old3 = old2
                old2 = old1
                old1 = work[i]
        #        print 'old1', old1
                if meth == 1:
                    #epsilon algorithm rule
                    work[i] = old3 + one / (work[i-1] - old2)
                else:
                    #quotient difference algorithm rules
                    if i % 2 == 0:
                        #difference form
       #                 print i, 'up', old3, work[i-1], old2
                        if old2 == old3:
      #                      print 'ouch', work[i-1]
                            work[i] = work[i-1]
                        elif old3 + work[i-1] == 0j:
     #                       print 'pow'
                            work[i] = old2
                        else:
                            work[i] = old3 - old2 + work[i-1] 
                    else:
                        #quotient form
    #                    print i, 'down', old3, work[i-1], old2
                        work[i] = old3 * work[i-1] / old2
   #                 print 'work[i-1] old2', work[i-1], old2, old3
                   # print 'work', work[i]
            #save continued fraction coefficients
            if meth != 1:
                d[j] = -1.0 * work[j]
   #             print 'd[j]', d[j]
 #   print 'W', work
  #  print 'd', d
    if meth == 1:
        #result of epsilon algorithm computation
        Result0 = work[N].real
    else:
        #evaluate continued fraction
        #initialize recurrence relations
        AOld1 = d[0]
        AOld2 = d[0]
        BOld1 = one + d[1] * Z
        BOld2 = one
        #use recurrence relations
        for j in range(2,N+1):
            if meth == 3 and j == N+1:
                #further acceleration on last iteration if required
                h2m = half * (one + (d[N] - d[N+1]) * Z)
                print (one + d[N+1] * Z / h2m**2)
                r2m = -1.0 * h2m * \
                      (one - sqrt((one + d[N+1] * Z / h2m**2)))
                A = AOld1 + r2m * AOld2
                B = BOld1 + r2m * BOld2
            else:
               # print 'd', j, d[j]
                A = AOld1 + d[j] * Z * AOld2
                #print 'A', A
                AOld2 = AOld1
                AOld1 = A
                B = BOld1 + d[j] * Z * BOld2
                #print 'B', B
                BOld2 = BOld1
                BOld1 = B
              #  print 'other', A, B
        ctemp = A /B
        #result of quotient difference algorithm evaluation
        Result0 = ctemp.real
    #return required approximate inverse
 #   print 'hhog', exp(gamma * t), Result0, bigT 
    return exp(gamma * t) * Result0 / bigT          
Back to Top

I am confused about how this code works as I am not an expert in python.
For example I do not know how I can define my function in Laplace domain in this code.

Lets say I need to find the inverse Laplace transform of the below function at t=1:

F(s) = 1/(s-1)

The inverse Laplace transform of F(s) is : e^(t)

At t=1 the result is expected to be = e

This code starts with :

def HoogTransform(t, gamma, bigT, N, FUN, meth, init, d, work):

I think the parameter FUN should be where I have to define the the F(s). I tried that but I could not get the code to work.

Could you please tell me how I can define F(s) in this code and how I can extract the result at t=1?
If it is not possible to get it to work and you know of any other simpler implementation of DE HOOG algorithm for inverse Laplace transform, I appreciate you to let me know.
Thanks in advance for your time and attention.

Best Regards
Jafar

RE: Laplase Inverse Transform Implementation - Larz60+ - Nov-15-2016

Please show your code with error

RE: Laplase Inverse Transform Implementation - jafar1363 - Nov-15-2016

This is the code with its documentation:

My Code

from numpy import *



def HoogTransform(t, gamma, bigT, N, FUN, meth, init, d, work):

'''Numerrical Laplace inversion using the Hoog method.

t (float variable) independent variable at which inverse

function is to be evaluated

gamma (float variable) parameter of inversion integral, governing

accuracy of inversion

bigT (float variable) parameter used to discretize inversion

integral, thus governing discretization error

N (integer variable) must be >= 2 and should be even; number

of terms to be used in sum for approximation, thus

governing truncation error

FUN name of user-supplied function which evaluates the complex

values of transform fbar(p) for complex values of Laplace

transform variable p

meth (integer variable) determines inversion method: 1=epsilon

method; 2=ordinary quotient difference algorithm;

3=acceleratedquotient difference algorithm

init (integer variable) =1 indicates this is the first call

of the algorithm with these values of gamma, bigT, N and

FUN. !=1 means array is not recomputed

d (complex array dimensioned in calling function) contains

continued fraction coefficient if init >1 and meth=2 or 3

work (complex array dimensioned in calling function) work space

used for calculating successive diagonals of the table



This is a translation of the FORTRAN subroutine LAPADC originally by

J.H. Knight. See code for more information.

Requires numpy.'''

#This is a translation of the FORTRAN subroutine LAPADC by J.H. Knight

#of the Division of Environmental Mechanics, CSIRO, Canberra, ACT2601,

#Australia originally coded in August 1986. The function performs a

#numerical LaPlace inversion using the Hoog method. Two papers are

#referenced in the original FORTRAN code:

#DeHoog, F.R., Knight, J.H. and Stokes, A.N. An improved method for

#numerical inversion of LaPlace transforms. SIAM J. SCI. STA. COMPUT.,

# 3, PP. 357-366, 1982.

#DeHoog, F.R., Knight, J.H. and Stokes, A.N. Subroutine LAPACC for LaPlace

#inversion by acceleration of convergence of complex sum. ACM TRANS.

#MATH. SOFTWARE, (manuscript in preperation).

#Translated by Karl W. Bandilla in July 2009.

#Comments are taken from the FORTRAN code



#Following is the header of the original FORTRAN code



# SUBROUTINE LAPADC(T,GAMMA,BIGT,N,F,FUN,METH,INIT,D,WORK)

#

#-----------------------------------------------------------------------

#

# AUTHOR

#

# J.H. KNIGHT

# DIVISION OF ENVIRONMENTAL MECHANICS

# CSIRO, CANBERRA, ACT 2601, AUSTRALIA

# AUGUST, 1986

#

# REFERENCES

#

# DE HOOG, F.R., KNIGHT, J.H. AND STOKES, A.N.

# AN IMPROVED METHOD FOR NUMERICAL INVERSION OF LAPLACE

# TRANSFORMS. SIAM J. SCI. STAT. COMPUT., 3, PP. 357-366, 1982.

#

# DE HOOG, F.R., KNIGHT, J.H. AND STOKES, A.N.

# SUBROUTINE LAPACC FOR LAPLACE TRANSFORM INVERSION BY

# ACCELERATION OF CONVERGENCE OF COMPLEX SUM.

# ACM TRANS. MATH. SOFTWARE, (MANUSCRIPT IN PREPARATION).

#

# ABSTRACT

#

# LAPADC COMPUTES AN APPROXIMATION TO THE INVERSE LAPLACE

# TRANSFORM F(T) CORRESPONDING TO THE LAPLACE TRANSFORM

# FBAR(P), USING A RATIONAL APPROXIMATION TO THE FOURIER

# SERIES RESULTING WHEN THE INVERSION INTEGRAL IS

# DISCRETIZED USING THE TRAPEZOIDAL RULE.

# THE RATIONAL APPROXIMATION IS A DIAGONAL PADE APPROXIMATION

# COMPUTED USING EITHER THE EPSILON ALGORITHM, IN WHICH THE

# EPSILON TABLE MUST BE RECOMPUTED AT EACH VALUE OF T, OR

# USING THE QUOTIENT-DIFFERENCE ALGORITHM, IN WHICH THE

# QUOTIENT-DIFFERENCE ALGORITHM IS USED TO COMPUTE THE

# COEFFICIENTS OF THE ASSOCIATED CONTINUED FRACTION, AND THE

# CONTINUED FRACTION IS EVALUATED AT EACH VALUE OF T.

# IN ADDITION, THE EVALUATION OF THE CONTINUED FRACTION IS

# MADE MORE ACCURATE WITH AN IMPROVED ESTIMATE OF THE

# REMAINDER.

#

# METHOD

#

# A REAL VALUED INVERSE FUNCTION F(T) IS THE REAL PART OF THE

# INTEGRAL

# INFINITY

#

# (ONE/PI)*EXP(GAMMA*T)* INTEGRAL FBAR(GAMMA+I*W)*EXP(I*W*T)*DW

#

# W=ZERO

#

# WITH FBAR THE LAPLACE TRANSFORM, GAMMA A REAL PARAMETER, AND

# I**2=-1. W IS THE VARIABLE OF INTEGRATION.

# IF THIS IS DISCRETIZED USING THE TRAPEZOIDAL RULE WITH STEP

# SIZE PI/BIGT, A FOURIER SERIES RESULTS, AND F(T) IS THE SUM

# OF THE DISCRETIZATION ERROR

#

# INFINITY

#

# SUM EXP(-TWO*GAMMA*K*BIGT)*F(T+TWO*K*BIGT)

#

# K=1

#

# AND THE REAL PART OF THE FOURIER SERIES

#

# INFINITY

#

# (ONE/BIGT)*EXP(GAMMA*T)* SUM A(K)*EXP(I*K*PI*T/BIGT)

#

# K=0

#

# WITH THE COEFFICIENTS (A(K),K=0,INFINITY) GIVEN BY

#

# A(0) = HALF*FBAR(GAMMA),

#

# A(K) = FBAR(GAMMA+I*K*PI/BIGT), K=1,INFINITY.

#

# THIS EXPRESSION CAN ALSO BE CONSIDERED AS A POWER SERIES IN

# THE VARIABLE Z=EXP(I*PI*T/BIGT), AND WE INTRODUCE A

# TRUNCATION ERROR WHEN WE USE A FINITE POWER SERIES

#

# M2

#

# V(Z,M2) = SUM A(K)*Z**K, M2 = M*2.

#

# K=0

#

# THE EPSILON ALGORITHM CALCULATES THE DIAGONAL PADE

# APPROXIMATION TO THIS FINITE POWER SERIES

#

# M M

#

# V(Z,M2) = SUM B(K)*Z**K / SUM C(K)*Z**K

#

# K=0 K=0

#

# + HIGHER ORDER TERMS, C(0)=ONE.

#

# THE CALCULATION MUST BE DONE FOR EACH NEW VALUE OF T.

# THIS METHOD IS USED FOR METH=1.

#

# THE QUOTIENT-DIFFERENCE (QD) ALGORITHM IS ANOTHER WAY

# OF CALCULATING THIS DIAGONAL PADE APPROXIMATION, IN THE FORM

# OF A CONTINUED FRACTION

#

# V(Z,M2) = D(0)/(ONE+D(1)*Z/(ONE+...+D(M2-1)*Z/(ONE+R2M(Z)))),

#

# WHERE R2M(Z) IS THE UNKNOWN REMAINDER.

# THE QD ALGORITHM GIVES THE COEFFICIENTS (D(K),K=0,M2), AND

# THE CONTINUED FRACTION IS THEN EVALUATED BY FORWARD RECURRENCE,

# FOR EACH NEW VALUE OF T, WITH R2M(Z) TAKEN AS D(M2)*Z.

# THE ANSWERS SHOULD BE IDENTICAL TO THOSE GIVEN BY THE

# EPSILON ALGORITHM, APART FROM ROUNDOFF ERROR.

# THIS METHOD IS USED FOR METH=2.

#

# IN ADDITION, THE CONVERGENCE OF THE CONTINUED FRACTION CAN

# BE ACCELERATED WITH AN IMPROVED ESTIMATE OF THE REMAINDER,

# SATISFYING

#

# R2M + ONE + (D(M2-1)-D(M2))*Z = D(M2)*Z/R2M

#

# THIS METHOD IS USED FOR METH=3.

#

# USAGE

#

# THE USER MUST CHOOSE THE PARAMETERS BIGT AND GAMMA WHICH

# DETERMINE THE DISCRETIZATION ERROR, AND N(=M2) WHICH

# DETERMINES THE TRUNCATION ERROR AND ROUNDOFF ERROR, GIVEN

# THE CHOICE OF BIGT AND GAMMA. AS N IS INCREASED, THE

# TRUNCATION ERROR DECREASES BUT THE ROUNDOFF ERROR INCREASES,

# SINCE THE EPSILON AND QD ALGORITHMS ARE NUMERICALLY UNSTABLE.

# IT IS DESIRABLE TO USE DOUBLE PRECISION FOR THE CALCULATIONS,

# UNLESS THE COMPUTER KEEPS AT LEAST 12 SIGNIFICANT DECIMAL

# DIGITS IN SINGLE PRECISION.

# METH=1 USES THE EPSILON ALGORITHM, METH=2 USES THE QUOTIENT-

# DIFFERENCE ALGORITHM, AND METH=3 USES THE QD ALGORITHM WITH

# IMPROVED ESTIMATE OF THE REMAINDER ON EVALUATION.

#

# AN ESTIMATE OF THE OPTIMAL VALUE OF N FOR GIVEN BIGT AND

# GAMMA CAN BE GOT BY VARYING N AND CHOOSING A VALUE SLIGHTLY

# LESS THAN THAT WHICH CAUSES THE RESULTS FOR METH=1 TO VARY

# SIGNIFICANTLY FROM THOSE FOR METH=2.

#

# CALCULATIONS USING THE QD ALGORITHM ARE MOST EFFICIENT

# AND MOST ACCURATE IF LAPADC IS CALLED WITH METH=3, AND WITH

# INIT=1 ON THE FIRST CALL WITH NEW VALUES OF BIGT, GAMMA, N

# OR FUN, AND WITH INIT=2 ON SUBSEQUENT CALLS WITH NEW VALUES

# OF T AND THE OTHER PARAMETERS UNCHANGED.

#

#

# DESCRIPTION OF ARGUMENTS

#

# T - REAL (DOUBLE PRECISION) VARIABLE. ON INPUT, VALUE

# OF INDEPENDENT VARIABLE AT WHICH INVERSE FUNCTION

# F(T) IS TO BE APPROXIMATED. UNCHANGED ON OUTPUT.

# GAMMA - REAL (DOUBLE PRECISION) VARIABLE. ON INPUT,

# CONTAINS PARAMETER OF INVERSION INTEGRAL, GOVERNING

# ACCURACY OF INVERSION. UNCHANGED ON OUTPUT.

# BIGT - REAL (DOUBLE PRECISION) VARIABLE. ON INPUT,

# CONTAINS PARAMETER USED TO DISCRETIZE INVERSION

# INTEGRAL. GOVERNS DISCRETIZATION ERROR.

# UNCHANGED ON OUTPUT.

# N - INTEGER VARIABLE, SHOULD BE EVEN .GE. 2. ON INPUT,

# CONTAINS NUMBER OF TERMS TO BE USED IN SUM FOR

# APPROXIMATION. DETERMINES TRUNCATION ERROR.

# UNCHANGED ON OUTPUT.

# F - REAL (DOUBLE PRECISION) VARIABLE. ON OUTPUT,

# CONTAINS VALUE OF INVERSE FUNCTION F(T), EVALUATED

# AT T.

# FUN - NAME OF USER-SUPPLIED SUBROUTINE WHICH EVALUATES

# (DOUBLE) COMPLEX VALUES OF TRANSFORM FBAR(P)

# FOR (DOUBLE) COMPLEX VALUES OF LAPLACE

# TRANSFORM VARIABLE P. CALLED IF METH .EQ. 1 .OR.

# INIT .EQ. 1. MUST BE DECLARED EXTERNAL IN

# (SUB)PROGRAM CALLING LAPADC.

# METH - INTEGER VARIABLE, DETERMINING METHOD TO BE USED TO

# ACCELERATE CONVERGENCE OF SUM.

# ON INPUT, METH .EQ. 1 INDICATES THAT THE EPSILON

# ALGORITHM IS TO BE USED.

# METH .EQ. 2 INDICATES THAT THE ORDINARY QUOTIENT

# DIFFERENCE ALGORITHM IS TO BE USED.

# METH .EQ. 3 INDICATES THAT THE QUOTIENT DIFFERENCE

# ALGORITHM IS TO BE USED, WITH FURTHER ACCELERATION

# OF THE CONVERGENCE OF THE CONTINUED FRACTION.

# UNCHANGED ON OUTPUT.

# INIT - INTEGER VARIABLE. ON INPUT, USED ONLY IF Q-D

# ALGORITHM IS SPECIFIED (METH .EQ. 2 OR 3).

# INIT .EQ. 1 INDICATES THAT THIS IS THE FIRST CALL

# OF THE ALGORITHM WITH THESE VALUES OF GAMMA, BIGT,

# N, FUN. IN THIS CASE, THE VALUES OF THE CONTINUED

# FRACTION COEFFICIENTS D(0:N) ARE RECALCULATED, AND

# THE CONTINUED FRACTION IS EVALUATED.

# INIT .NE. 1 INDICATES THAT THE VALUES OF GAMMA,

# BIGT, N, AND FUN ARE UNCHANGED SINCE THE LAST CALL.

# IN THIS CASE, THE INPUT VALUES OF THE CONTINUED

# FRACTION COEFFICIENTS D(0:N) ARE USED, AND THE

# CONTINUED FRACTION IS EVALUATED.

# UNCHANGED ON OUTPUT.

# D - (DOUBLE) COMPLEX ARRAY, OF VARIABLE DIMENSION (0:N).

# MUST BE DIMENSIONED IN CALLING PROGRAM.

# ON INPUT, IF METH .EQ. 2 OR 3 .AND. INIT .NE. 1,

# MUST CONTAIN CONTINUED FRACTION COEFFICIENTS

# OUTPUT ON A PREVIOUS CALL.

# ON OUTPUT, IF METH .EQ. 2 OR 3, CONTAINS CONTINUED

# FRACTION COEFFICIENTS CALCULATED ON THIS CALL

# (INIT .EQ. 1) OR AS SUPPLIED ON INPUT (INIT .NE. 1).

# WORK - (DOUBLE) COMPLEX ARRAY OF VARIABLE DIMENSION (0:N).

# MUST BE DIMENSIONED IN CALLING PROGRAM.

# WORK SPACE USED FOR CALCULATING SUCCESSIVE DIAGONALS

# OF THE TABLE IF METH .EQ. 1 .OR. INIT .EQ. 1.

# ON OUTPUT, CONTAINS LAST CALCULATED DIAGONAL OF

# TABLE.

#

# SUBROUTINE CALLED

#

# FUN(P,FBAR) - USER-SUPPLIED SUBROUTINE WHICH EVALUATES

# (DOUBLE) COMPLEX VALUES OF TRANSFORM FBAR(P)

# FOR (DOUBLE) COMPLEX VALUES OF LAPLACE

# TRANSFORM VARIABLE P. CALLED IF METH .EQ. 1 .OR.

# INIT .EQ. 1. MUST BE DECLARED EXTERNAL IN

# (SUB)PROGRAM WHICH CALLS SUBROUTINE LAPADC.





if N < 2:

print 'Order too low (N at least 2)'

return -1.0



zeror = 0.0

zero = 0.0 + 0.0j

half = 0.5 + 0.0j

one = 1.0 + 0.0j

# work = zeros(N+1,dtype=complex128)

# d = zeros(N+1,dtype=complex128)

#make order multiple of 2

M2 = (N / 2) * 2

factor = pi / bigT

# print 'factor', factor

argR = t * factor

Z = cos(argR) + sin(argR) * 1j

if meth == 1 or init == 1:

#new: get all A's first

##########################

#need to think of way to do this

###########################

#calculate table if epsilon algorithm is to be used, or

#if this is the first call to quotient-difference algorithm

#with these parameters and inverse transform

A = FUN(gamma + zeror * 1j)

AOld = half * A

if AOld == 0.0:

return 0.0

argI = factor

A = FUN(gamma + argI * 1j)

#initialize table entries

if meth == 1:

ztoj = Z

term = A * ztoj

work[0] = AOld + term

work[1] = one / term

else:

d[0] = AOld

work[0] = zero

work[1] = A / AOld

d[1] = -1.0 * work[1]

AOld = A

#calculate successive diagonals of table

for j in range(2, N + 1):

#initialize calculation of diagonal

old2 = work[0]

old1 = work[1]

argI += factor

A = FUN(gamma + argI * 1j)

# print 'A', A, A - AOld

if A == 0:

return 0.0

if meth == 1:

#calculate next term and sum of power series

ztoj = Z * ztoj

term = A * ztoj

work[0] += term

work[1] = one / term

else:

work[0] = zero

work[1] = A / AOld

AOld = A

# print work[0:10]

#calculate diagonal using Rhombus rules

for i in range(2, j+1):

old3 = old2

old2 = old1

old1 = work[i]

# print 'old1', old1

if meth == 1:

#epsilon algorithm rule

work[i] = old3 + one / (work[i-1] - old2)

else:

#quotient difference algorithm rules

if i % 2 == 0:

#difference form

# print i, 'up', old3, work[i-1], old2

if old2 == old3:

# print 'ouch', work[i-1]

work[i] = work[i-1]

elif old3 + work[i-1] == 0j:

# print 'pow'

work[i] = old2

else:

work[i] = old3 - old2 + work[i-1] 

else:

#quotient form

# print i, 'down', old3, work[i-1], old2

work[i] = old3 * work[i-1] / old2

# print 'work[i-1] old2', work[i-1], old2, old3

# print 'work', work[i]

#save continued fraction coefficients

if meth != 1:

d[j] = -1.0 * work[j]

# print 'd[j]', d[j]

# print 'W', work

# print 'd', d

if meth == 1:

#result of epsilon algorithm computation

Result0 = work[N].real

else:

#evaluate continued fraction

#initialize recurrence relations

AOld1 = d[0]

AOld2 = d[0]

BOld1 = one + d[1] * Z

BOld2 = one

#use recurrence relations

for j in range(2,N+1):

if meth == 3 and j == N+1:

#further acceleration on last iteration if required

h2m = half * (one + (d[N] - d[N+1]) * Z)

print (one + d[N+1] * Z / h2m**2)

r2m = -1.0 * h2m * \

(one - sqrt((one + d[N+1] * Z / h2m**2)))

A = AOld1 + r2m * AOld2

B = BOld1 + r2m * BOld2

else:

# print 'd', j, d[j]

A = AOld1 + d[j] * Z * AOld2

#print 'A', A

AOld2 = AOld1

AOld1 = A

B = BOld1 + d[j] * Z * BOld2

#print 'B', B

BOld2 = BOld1

BOld1 = B

# print 'other', A, B

ctemp = A /B

#result of quotient difference algorithm evaluation

Result0 = ctemp.real

#return required approximate inverse

# print 'hhog', exp(gamma * t), Result0, bigT 

return exp(gamma * t) * Result0 / bigT

The thing is I have to define my function in the Laplace domain somewhere in this code. For example:

F(s) = 1/(s-1)

This code computes the inverse of F(s) which is :

f(t)=e^t = e (if t=1)

I don't know how I can define F(s) and extract the results for example at t=1.
Thanks in advance for your help.

Best
Jafar

RE: Laplase Inverse Transform Implementation - Larz60+ - Nov-15-2016

Ouch! Much more code than wanted. Is it possible to make a small snippet that will run and reproduce the error?
Looks like all indentation was lost
This shouldn't happer if you cut and paste fro, a properly formatted python script if you use cut and paste.

RE: Laplase Inverse Transform Implementation - sparkz_alot - Nov-15-2016

when pasting, use [control] + [shift] + 'v'

yes, really don't need to see ~600 lines of comments. Please post the entire Traceback error as well.

RE: Laplase Inverse Transform Implementation - jafar1363 - Nov-15-2016

Thanks for your responses.
I do not know how familiar you are with Laplace inverse transform. First of all we need to define a function in Laplace domain like:

F(s) = 1 / (s-1)

After implementation of Laplace inverse we will have a function in the time-domain. For example the Laplace inverse of the above function is :

f(t) = e^t

I have seen a lot of other python codes (based on other algorithms) and in all of them there is a certain part for defining the F(s). For example:

#   The Laplace function F(s) is defined here.
def F(s):
    return 1.0 / (s - 1.0)

But there is no part in this code that the user can define the F(s). That is why I can not even run this code. (How am I supposed to get a error when I cant run the code?!).
Could you please take a closer look at this code and figure out how we can get it to work? Thanks again for your time.

Jafar

RE: Laplase Inverse Transform Implementation - jafar1363 - Nov-15-2016

I recently added some lines at the bottom of that code to define the F(s) and call the HoogTransform function from the code with its parameters. These are the 3 new lines I added at the end of the code:

#defining the F(s)=1/(s-1)
def FUN(s):
   return 1.0/(s-1.0)

#Calling the HoogTransform function
HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10)

And this is the complete code:

Hide/Show

from numpy import *



def HoogTransform(t, gamma, bigT, N, FUN, meth, init, d, work):

   if N < 2:
       print 'Order too low (N at least 2)'
       return -1.0

   zeror = 0.0
   zero = 0.0 + 0.0j
   half = 0.5 + 0.0j
   one = 1.0 + 0.0j
#    work = zeros(N+1,dtype=complex128)
#    d = zeros(N+1,dtype=complex128)
   #make order multiple of 2
   M2 = (N / 2) * 2
   factor = pi / bigT
#   print 'factor', factor
   argR = t * factor
   Z = cos(argR) + sin(argR) * 1j
   if meth == 1 or init == 1:
       #new: get all A's first
       ##########################
       #need to think of way to do this
       ###########################
       #calculate table if epsilon algorithm is to be used, or
       #if this is the first call to quotient-difference algorithm
       #with these parameters and inverse transform
       A = FUN(gamma + zeror * 1j)
       AOld = half * A
       if AOld == 0.0:
           return 0.0
       argI = factor
       A = FUN(gamma + argI * 1j)
       #initialize table entries
       if meth == 1:
           ztoj = Z
           term = A * ztoj
           work[0] = AOld + term
           work[1] = one / term
       else:
           d[0] = AOld
           work[0] = zero
           work[1] = A / AOld
           d[1] = -1.0 * work[1]
           AOld = A
       #calculate successive diagonals of table
       for j in range(2, N + 1):
           #initialize calculation of diagonal
           old2 = work[0]
           old1 = work[1]
           argI += factor
           A = FUN(gamma + argI * 1j)
        #   print 'A', A, A - AOld
           if A == 0:
               return 0.0
           if meth == 1:
               #calculate next term and sum of power series
               ztoj = Z * ztoj
               term = A * ztoj
               work[0] += term
               work[1] = one / term
           else:
               work[0] = zero
               work[1] = A / AOld
               AOld = A
       #    print work[0:10]
           #calculate diagonal using Rhombus rules
           for i in range(2, j+1):
               old3 = old2
               old2 = old1
               old1 = work[i]
       #        print 'old1', old1
               if meth == 1:
                   #epsilon algorithm rule
                   work[i] = old3 + one / (work[i-1] - old2)
               else:
                   #quotient difference algorithm rules
                   if i % 2 == 0:
                       #difference form
      #                 print i, 'up', old3, work[i-1], old2
                       if old2 == old3:
     #                      print 'ouch', work[i-1]
                           work[i] = work[i-1]
                       elif old3 + work[i-1] == 0j:
    #                       print 'pow'
                           work[i] = old2
                       else:
                           work[i] = old3 - old2 + work[i-1]
                   else:
                       #quotient form
   #                    print i, 'down', old3, work[i-1], old2
                       work[i] = old3 * work[i-1] / old2
  #                 print 'work[i-1] old2', work[i-1], old2, old3
                  # print 'work', work[i]
           #save continued fraction coefficients
           if meth != 1:
               d[j] = -1.0 * work[j]
  #             print 'd[j]', d[j]
#   print 'W', work
 #  print 'd', d
   if meth == 1:
       #result of epsilon algorithm computation
       Result0 = work[N].real
   else:
       #evaluate continued fraction
       #initialize recurrence relations
       AOld1 = d[0]
       AOld2 = d[0]
       BOld1 = one + d[1] * Z
       BOld2 = one
       #use recurrence relations
       for j in range(2,N+1):
           if meth == 3 and j == N+1:
               #further acceleration on last iteration if required
               h2m = half * (one + (d[N] - d[N+1]) * Z)
               print (one + d[N+1] * Z / h2m**2)
               r2m = -1.0 * h2m * \
                     (one - sqrt((one + d[N+1] * Z / h2m**2)))
               A = AOld1 + r2m * AOld2
               B = BOld1 + r2m * BOld2
           else:
              # print 'd', j, d[j]
               A = AOld1 + d[j] * Z * AOld2
               #print 'A', A
               AOld2 = AOld1
               AOld1 = A
               B = BOld1 + d[j] * Z * BOld2
               #print 'B', B
               BOld2 = BOld1
               BOld1 = B
             #  print 'other', A, B
       ctemp = A /B
       #result of quotient difference algorithm evaluation
       Result0 = ctemp.real
   #return required approximate inverse
#   print 'hhog', exp(gamma * t), Result0, bigT
   return exp(gamma * t) * Result0 / bigT

def FUN(s):
   return 1.0/(s-1.0)

HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10)

But I get this error:

Error:Traceback (most recent call last):
  File "/home/hp/Desktop/4.py", line 147, in <module>
    HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10)
  File "/home/hp/Desktop/4.py", line 46, in HoogTransform
    d[0] = AOld
TypeError: 'int' object does not support item assignment

Does anybody have any idea about how I can fix it and finally get the code to work?
Thanks

RE: Laplase Inverse Transform Implementation - Larz60+ - Nov-16-2016

HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10) ???
Where is this function?

RE: Laplase Inverse Transform Implementation - jafar1363 - Nov-16-2016

(Nov-16-2016, 12:07 AM)Larz60+ Wrote: HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10) ???
Where is this function?

Just click on the spoiler in my previous reply. Go to the line 5. This the function that returns the Laplace inverse of FUN. As you can see it includes different parameters. I just substitute the parameters with some arbitrary values. In addition I defined the FUN but it did not work. I am still confused about how this code works!

RE: Laplase Inverse Transform Implementation - nilamo - Nov-17-2016

(Nov-15-2016, 11:55 PM)jafar1363 Wrote: But I get this error:
Error:Traceback (most recent call last):
  File "/home/hp/Desktop/4.py", line 147, in <module>
    HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10)
  File "/home/hp/Desktop/4.py", line 46, in HoogTransform
    d[0] = AOld
TypeError: 'int' object does not support item assignment
Does anybody have any idea about how I can fix it and finally get the code to work?
Thanks

>>> d = 42
>>> d[0] = "spam"
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'int' object does not support item assignment

If you want "d" to be a list, you should make it a list, instead of the int it actually is.