Laplace Inverse Transform Implementation

[jafar1363]

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