井流模型的Laplace变换数值反演
有了井流模型的Laplace 变换解后，大部分情况下很难通过 Laplace 逆变换得出解的解析形式，只有在有限的情况可以容易得出解析。这时 Laplace 变换数值反演就有了用武之地。
以下内容源自地下水动力学课程的补充材料，以 Theis 解为例，其他井流函数可以用类似方法求得。
式中，。
取 ，得
由 Laplace 变换性质
记 ，
，有
取 ，即
，可以计算出
。
Stehfest 算法
Stehfest 方法是一种非常简单的拉普拉斯逆变换数值方法。该算法速度快，通常能给出很好的结果，尤其是对于光滑函数。值得注意的是，Stehfest 首次发表的数学模型是错误的，后来修订如下：
设 为函数
的 Laplace 变换，Stehfest 算法公式如下：
式中， 为 Stehfest 系数，按如下公式计算：
先定义一个测试函数：
def timer(func):
def func_wrapper(*args, **kwargs):
from time import time
time_start = time()
result = func(*args, **kwargs)
time_end = time()
time_spend = time_end - time_start
print("%s cost time: %.6f s" % (func.__name__, time_spend))
return result
return func_wrapper
Stehfest 算法程序
import math as math
import scipy.special as sp
'''
This uses the Stehfest algorithm to invert a Laplace transform. See Stehfest,H.
(1970) "Numerical inversion of Laplace transforms," Comm. ACM, Vol.13, No.1,
pages 47-49 and 624. The approximation has the nature of an asymptotic series,
in which n is the number of terms used. Thus, n shuld not exceed about 20, and
n = 10 is probably sufficient for most applications. Note that n must be an even
integer and that t = time. Beware: this works well for some transforms but very
poorly for others. Furthermore, a value for n that is slightly too large can
cause a dramatic decrease in accuracy.
Copyright © Dr Bruce Hunt 2003
'''
def stehcoef(i, n):
nhlf = int(n/2)
k1 = (i + 1)//2
k2 = min(i, nhlf)
stehcoef = 0.0
for k in range(k1, k2 + 1):
stehcoef = stehcoef + k**nhlf*math.factorial(2*k)/math.factorial(nhlf-k)/ \
math.factorial(k)/math.factorial(k-1)/math.factorial(i-k)/math.factorial(2*k-i)
return (-1)**(i + nhlf)*stehcoef
def lapinv(t, N):
if N%2 != 0: N += 1
lapinv = 0.0
a = math.log(2.0)/t
for i in range(1,N+1):
lapinv = lapinv+stehcoef(i, N)*transform(i*a)
return lapinv*a
验证：
def transform(p):
return 2*sp.kv(0, 2*math.sqrt(p))/p
@timer
def test_Stehfest():
u = []
W = []
for i in range(-5, 2):
u.append(10**i)
W.append(lapinv(10**(-i), 20))
print(u)
print(W)
test_Stehfest()
[1e-05, 0.0001, 0.001, 0.01, 0.1, 1, 10]
[10.935586387181772, 8.633328189684693, 6.331430232420058, 4.037790209571647, 1.8229829978456118, 0.21938416309114728, 4.189303016734773e-06]
test_Stehfest cost time: 0.000685 s
Talbot 算法
刚开始学习时看的是英文文献，意思是逆变换等价于复平面 Bromwich 轮廓积分，如果轮廓为围绕负实轴的抛物线，Talbot 方法可使积分快速收敛。想进一步了解 Talbot 算法可参考相关文献。
Talbot 算法用到复数计算库 cmath中的 cmath.exp()函数。
Talbot 算法程序
# -*- coding: utf-8 -*-
'''
Talbot suggested that the Bromwich line be deformed into a contour that begins
and ends in the left half plane, i.e., z → −∞ at both ends.
Due to the exponential factor the integrand decays rapidly
on such a contour. In such situations the trapezoidal rule converge
extraordinarily rapidly.
For example here we compute the inverse transform of F(s) = 1/(s+1) at t = 1
>>> error = Talbot(1,24)-exp(-1)
>>> error
(3.3306690738754696e-015+0j)
Talbot method is very powerful here we see an error of 3.3e-015
with only 24 function evaluations
Created by Fernando Damian Nieuwveldt
email:fdnieuwveldt@gmail.com
Date : 25 October 2009
Reference
L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
and rational approximations. BIT. Numerical Mathematics,
46(3):653 670, 2006.
'''
import scipy.special as sp # scipy.special 库计算Bessel函数，如 sp.kv(n, x)等。
import cmath as cmath # 程序中用 cmath 库实现复数运算。
def talbot(t, N):
# Initiate the stepsize
h = 2.0*cmath.pi/N;
# Shift contour to the right in case there is a pole on the positive real axis :
# Note the contour will not be optimal since it was originally devoloped for function
# with singularities on the negative real axis
# For example take F(s) = 1/(s-1), it has a pole at s = 1, the contour needs to be shifted
# with one unit, i.e shift = 1. But in the test example no shifting is necessary
shift = 0.0;
ans = 0.0;
if t == 0:
print("ERROR: Inverse transform can not be calculated for t=0")
return ("Error")
# The for loop is evaluating the Laplace inversion at each point theta which is based on the
# trapezoidal rule
for k in range(0, N):
theta = -cmath.pi + (k + 0.5)*h
z = shift + N/t*(0.5017*theta/cmath.tan(0.6407*theta) - 0.6122 + 0.2645j*theta)
dz = N/t*(-0.5017*0.6407*theta/(cmath.sin(0.6407*theta)**2) + 0.5017/cmath.tan(0.6407*theta) + 0.2645j)
ans = ans + cmath.exp(z*t)*transform(z)*dz;
return ((h/(2.0j*cmath.pi))*ans).real
验证：
# Here is the Laplace function to be inverted, should be changed manually
def transform(p):
return 2*sp.kv(0, 2*cmath.sqrt(p))/p
N = 24
@timer
def test_talbot():
u = []
W = []
for i in range(-5, 2):
u.append(10**i)
W.append(talbot(10**(-i), N))
print(u)
print(W)
test_talbot()
[1e-05, 0.0001, 0.001, 0.01, 0.1, 1, 10]
[10.935719800043493, 8.63322470457448, 6.331539364135962, 4.037929576537966, 1.8229239584192654, 0.21938393439543072, 4.156968879129344e-06]
test_talbot cost time: 0.001674 s
用于任意精度浮点运算的 Python 库 mpmath 中也有 Laplace 数值反演的函数 mpmath.invertlaplace()，需要注意的是程序中所用的 Bessel 函数也必须用库中的函数 mpmath.besselk()。
我们知道，Theis 解的 Laplace 变换为如下形式
式中，。
井函数还可以按如下方法计算。
- 按上式定义 Laplace 变换函数；
- 给出参数
的值，如都取值 1；
- 给出
的值，如
u=[1e-05, 1e-04, 1e-03, 1e-02, 1e-01, 1, 10]； - 按
计算
；
- 调用
talbot(t, N)函数，返回值为降深；
- 井函数的值为
。
代码如下：
Q, S, T, a, r = 1., 1., 1., 1., 1.
def transform(p):
return Q*sp.kv(0, cmath.sqrt(p/a)*r)/(2.0*math.pi*T*p)
N = 24
@timer
def test_talbot2():
u = []
W = []
for i in range(-5, 2):
u1 = 10**i
u.append(u1)
t = r**2.*S/(4.*T*u1)
W.append(talbot(t, N)*4.*math.pi*T/Q)
print(u)
print(W)
test_talbot2()
[1e-05, 0.0001, 0.001, 0.01, 0.1, 1, 10]
[10.935719800043474, 8.63322470457447, 6.331539364135964, 4.037929576537967, 1.8229239584192658, 0.21938393439543072, 4.156968879129344e-06]
test_talbot2 cost time: 0.001005 s
对比计算结果可以得出，Talbot 算法比 Stehfest 算法精度高，对比其他方法数值反演更为灵活。
