井流模型的Laplace變換數(shù)值反演
有了井流模型的Laplace 變換解后，大部分情況下很難通過 Laplace 逆變換得出解的解析形式毒返，只有在有限的情況可以容易得出解析。這時 Laplace 變換數(shù)值反演就有了用武之地。
以下內(nèi)容源自地下水動力學(xué)課程的補(bǔ)充材料凡辱，以 Theis 解為例，其他井流函數(shù)可以用類似方法求得栗恩。
式中透乾，。
取 ，得
由 Laplace 變換性質(zhì)
記 再来，
，有
取 磷瘤，即
芒篷，可以計算出
。
Stehfest 算法
Stehfest 方法是一種非常簡單的拉普拉斯逆變換數(shù)值方法膀斋。該算法速度快梭伐，通常能給出很好的結(jié)果，尤其是對于光滑函數(shù)仰担。值得注意的是糊识，Stehfest 首次發(fā)表的數(shù)學(xué)模型是錯誤的，后來修訂如下：
設(shè) 為函數(shù)
的 Laplace 變換摔蓝，Stehfest 算法公式如下：
式中赂苗， 為 Stehfest 系數(shù)，按如下公式計算：
先定義一個測試函數(shù)：
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 算法
剛開始學(xué)習(xí)時看的是英文文獻(xiàn)贮尉，意思是逆變換等價于復(fù)平面 Bromwich 輪廓積分拌滋，如果輪廓為圍繞負(fù)實軸的拋物線，Talbot 方法可使積分快速收斂猜谚。想進(jìn)一步了解 Talbot 算法可參考相關(guān)文獻(xiàn)败砂。
Talbot 算法用到復(fù)數(shù)計算庫 cmath中的 cmath.exp()函數(shù)。
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函數(shù)魏铅，如 sp.kv(n, x)等昌犹。
import cmath as cmath # 程序中用 cmath 庫實現(xiàn)復(fù)數(shù)運算。
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 數(shù)值反演的函數(shù) mpmath.invertlaplace()览芳，需要注意的是程序中所用的 Bessel 函數(shù)也必須用庫中的函數(shù) mpmath.besselk()斜姥。
我們知道，Theis 解的 Laplace 變換為如下形式
式中沧竟，铸敏。
井函數(shù)還可以按如下方法計算。
- 按上式定義 Laplace 變換函數(shù)悟泵；
- 給出參數(shù)
的值杈笔，如都取值 1；
- 給出
的值糕非，如
u=[1e-05, 1e-04, 1e-03, 1e-02, 1e-01, 1, 10]蒙具； - 按
計算
敦第；
- 調(diào)用
talbot(t, N)函數(shù)，返回值為降深店量；
- 井函數(shù)的值為
芜果。
代碼如下：
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
對比計算結(jié)果可以得出，Talbot 算法比 Stehfest 算法精度高融师，對比其他方法數(shù)值反演更為靈活右钾。
