Posts: 6
Threads: 1
Joined: Nov 2016
Nov-15-2016, 03:59 AM
(This post was last modified: Nov-17-2016, 09:03 PM by Ofnuts.)
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:
#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
Posts: 11,873
Threads: 474
Joined: Sep 2016
Please show your code with error
Posts: 6
Threads: 1
Joined: Nov 2016
Nov-15-2016, 05:36 PM
(This post was last modified: Nov-15-2016, 06:00 PM by metulburr.)
This is the code with its documentation:
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
Posts: 11,873
Threads: 474
Joined: Sep 2016
Nov-15-2016, 05:51 PM
(This post was last modified: Nov-15-2016, 05:52 PM by Larz60+.)
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.
Posts: 1,298
Threads: 38
Joined: Sep 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.
If it ain't broke, I just haven't gotten to it yet.
OS: Windows 10, openSuse 42.3, freeBSD 11, Raspian "Stretch"
Python 3.6.5, IDE: PyCharm 2018 Community Edition
Posts: 6
Threads: 1
Joined: Nov 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
Posts: 6
Threads: 1
Joined: Nov 2016
Nov-15-2016, 11:55 PM
(This post was last modified: Nov-16-2016, 12:14 AM by metulburr.)
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:
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
Posts: 11,873
Threads: 474
Joined: Sep 2016
HoogTransform (1.0, 100, 20, 100, FUN, 3, 1, 10, 10) ???
Where is this function?
Posts: 6
Threads: 1
Joined: Nov 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!
Posts: 3,458
Threads: 101
Joined: Sep 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.
|
