Page 12 - Laplace transform numerical inversion

Bruno Josso & Leif Larsen: Laplace transform numerical inversion - June 2012 - p 10/18
2 3 4 5 6
-0.5
0 1
0
0 1 2 3 4 5 6
-0.003
-0.002
-0.001
0
0.001
0.002
0.003
Error: f(t) - Fp_inv(t)
Test: Fp_inv(t) < 0
Figure 8: Inversion of f (t) = e− µt by Stehfest N = 22 algorithm.
In a standard implementation, the Stehfest number N should be adjustable. One can notice here
that it is reasonable to only allow a Stehfest number N ≤ 18. Indeed, if one tries to run Ste-
hfest algorithm with N bigger than 18 by using standard type formats (up to "long double" real
data type), one will suffer numerical approximation errors. Hence, when the Stehfest number N
becomes too big, the result is opposite to what was expected: reconstruction becomes less accu-
rate. To overcome this numerical approximation problem, one has to change the implementation
method by using: "arbitrary precision floating point arithmetic". With such an implementation
Gaver-Stehfest algorithm can then be pushed further. It can then be called: "big number Gaver-
Stehfest algorithm".
4.1.2 Big number Gaver-Stehfest algorithm
What we call the "big number Gaver-Stehfest algorithm" corresponds to the Gaver-Stehfest in-
version algorithm for which the Stehfest number N typically becomes greater than 18. One could
see in section( 4.1.1) that though in standard implementation "long double" type is coded with
80 bits, numerical approximation errors occur starting from N = 18. Now, if one wants to go
further with running the algorithm with N > 18, one needs to use arbitrary precision floating point
arithmetic.
Even when changing the floating point arithmetic, Gaver-Stehfest algorithm is implemented ex-
actly the same way as before. One keeps using the approximation(18) to inverse the Laplace
transform. It should also be pointed out that in the reservoir engineering context, one also needs
to reimplement the reservoir models by also using the arbitrary precision floating point arithmetic.
Only changing the algorithm but keeping the reservoir models old implementation does not solve
the numerical approximation problem.
One can see in Figure 9,the big number Gaver-Stehfest algorithm reconstruction errors and sign
test for N = 22 , 26, 30 and 40 (to be compared with graphes from Figure 4 to Figure 8).
One can see in Figure 9 that with such an implementation, reconstruction behaviour becomes as
expected: the bigger N the more accurate the reconstruction and the smaller the error. It can
virtually be made as small as desired. It can also be noticed that reconstruction error is smooth;
direct comparison can be made for N = 22. This means that, with such an arbitrary precision
floating point arithmetic implementation, results are not impacted by numerical precision any
more.
When checking sign test curves in green in Figure 9,one can see that even if reconstruction error
decreases with N , one can never avoid the inverted function to become negative. Depending on
N , the negative swing occurs at different instants. Unfortunately, this is not straightforwardly
controllable and therefore cannot easily be rejected towards very late time. Generalising this
result, one can say that for all retained Stehfest number N , the reconstructed function always