拐点法的程序（Python）
欢迎点评！ 如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。
需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions
用 pip 安装：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
以下内容源自地下水动力学课程教学内容，计算井函数的程序要独立保存为 wellfunction.py，程序可在 Jupyter notebook 中运行。
拐点法是根据 Hantush-Jacob 公式得出的 曲线拐点的特殊性质求参数。
根据 Hantush-Jacob 公式
拐点 的性质
- 坐标与降深：
- 切线斜率：
-
、
的关系：
思路
- 调整参数确定拐点坐标与降深
，最大降深
用外推法确定;
- 调整参数确定拐点切线斜率
- 二分法求解方程:
:
是单调递减的，二分法区间端点函数值符号应相反，使二分法有效。
- 计算参数：
- 验证：避免确定
程序先根据原始数据确定最大降深 sp，并估计出 tp，用 tp 相邻的点计算初始的 T,S,slope，通过调节 tp,sp,slope 获得水文地质参数。
上代码：
%matplotlib widget
import numpy as np
import math as math
from scipy.optimize import bisect
from scipy.special import k0
import matplotlib.pyplot as plt
import ipywidgets as widgets
from wellfunction import *
# 控制小数的显示精度
np.set_printoptions(precision=4)
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
# 准备数据
Q = 69.1/60 # m^3/min
r = 197.0 # m
# min
t = np.array([
1, 4, 7, 10, 15, 20, 25, 30, 45, 60, \
75, 90, 120, 150, 180, 210, 240, 270, 300, 330, \
360, 390, 420, 450, 480, 510, 540])
# m
s = np.array([
0.05, 0.054, 0.12, 0.175, 0.26, 0.33, 0.383, 0.425, 0.52, \
0.575, 0.62, 0.64, 0.685, 0.725, 0.735, 0.755, 0.76, 0.76, \
0.763, 0.77, 0.772, 0.785, 0.79, 0.792, 0.794, 0.795,
0.796])
n = len(t)
# 计算绘图范围
imin = math.floor(math.log10(min(t))) # math.floor(x)返回小于x的最大整数
imax = math.ceil(math.log10(max(t))) # math.ceil(x)返回大于等于x的最小整数
xmin = 10**imin
xmax = 10**imax
ymin = 0.0
ymax = math.ceil(max(s*10))/10
# 绘函数图像的网格
x = np.linspace(imin, imax, (imax-imin)*10+1)
x = np.float_power(10,x)
# 设定初始的 sp, tp
sp = 0.5*s[-1]
idx = np.where(s > sp)[0][0]
tp = t[idx]
# 设置初始的 T, S, beta, i
# 初始参数，最简单方法是取 tp 相邻点用 Jacob 公式计算
n = len(t)
i1 = idx
i2 = i1 + 1
t1 = t[i1]
t2 = t[i2]
s1 = s[i1]
s2 = s[i2]
kk = (s1 - s2)/np.log10(t1/t2)
T = 0.183*Q/kk
S = 2.25*T*t1/r**2/np.float_power(10, s1/kk)
B = 999.0
slope = kk
# Plot the data
def inflection_point_fit(tp, sp, slope):
global T, S, B # 这些都是全局变量，改变后函数外的值同时改变
plt.style.use('default') # 绘图风格
fig = plt.figure()
ax = fig.add_subplot()
'''
若不想设置其他内容，也可合并写为
fig, ax = plt.subplots(figsize=(6, 4))
'''
ax.plot(t, s, '*', label="观测值")
ax.plot(x, sp + slope*np.log10(x/tp), label="拐点切线")
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_xscale("log")
plt.xlabel('$\log t$')
plt.ylabel('$s$')
ax.grid(True)
ax.set_title("拐点法", fontproperties={'family': 'KaiTi', 'size': 12}) # 指定显示中文字体
# 调用 scipy.optimize.bisect(二分法)
beta = bisect(lambda x:np.exp(x)*k0(x)-2.3*sp/slope,0.01,5)
# 计算参数
T = 2.3*Q*np.exp(-beta)/4/np.pi/slope
S = 2*T*tp*beta/r**2
B = r/beta
# 绘图数据
u = 0.25*r**2*S/T/x
ax.plot(x, 0.25*Q*hantush_jacob(u, beta)/np.pi/T, label="标准曲线")
plt.legend(
prop={'family': 'Simsun', 'size': 10}, handlelength=2,
loc=4,title="图例",
title_fontproperties={'family': 'KaiTi', 'size': 12})
plt.show()
print(' T(m^2/min): ', '{:.4f}'.format(T))
print(' S: ', '{:.4e}'.format(S))
print(' B(m): ', '{:.4e}'.format(B))
widgets.interact(
inflection_point_fit,
tp = widgets.FloatSlider(
value=tp, min=1, max=0.75*t[-1], step=1,
description=r'$t_p$ [-]:', continuous_update=False,
readout_format='.1f', disabled=False),
sp = widgets.FloatText(
value=sp, description=r'$s_p$ [-]:',
continuous_update=False, disabled=False),
slope = widgets.FloatSlider(
value=slope, min=0.1, max=1, step=0.01,
description=r'$slope$ [-]:', continuous_update=False,
readout_format='.3f', disabled=False)
);
print(' T(m^2/min): ', '{:.4f}'.format(T))
print(' S: ', '{:.4e}'.format(S))
print(' B(m): ', '{:.4e}'.format(B))
image.png
