配线法的程序实现（Python）
欢迎点评！如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。
需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions
用 pip 安装：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
ipywidgets 模块可实现交互绘图。
以下内容源自地下水动力学课程教学内容，计算井函数的程序要独立保存为 wellfunction.py，程序可在 Jupyter notebook 中运行。
配线法原理
Theis 配线法是 Hantush-Jacob 配线法的特例，在此只介绍后一种。
对
两边取对数
思路：
在标准曲线坐标系中移动 散点图，根据平移量计算
：
记坐标平移量为
通过输入确定。
参数计算公式：
代码中给出一些关于初始参数的解释。
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
import math as math
from wellfunction import *
# 控制小数的显示精度
np.set_printoptions(precision=4)
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
# 准备数据, 根据配线法的种类，需要专门准备数据
Q = 528/1440 # m^3/min
r = 90 # m
t = np.array([1, 2, 4, 6, 9, 20, 30, 40, 50, 60, 90, 120, 150, 360, 550, 720]) # min
s = np.array([2.5, 3.9, 6.1, 8.0, 10.6, 16.8, 20.0, 22.6, 24.7, 26.4, 30.4, 33.0, 35.0, 42.6, 44.0, 44.5])/100 # m
ymin = 1.0e-2
ymax = 10.0
xmin = 1.0e-1
xmax = 1.0e4
x = np.linspace(-1, 4, 51)
x = np.float_power(10,x)
# 初始参数，最简单方法是取中间相邻两点用 Jacob 公式计算
# 取数组长度
n = len(t)
i1 = int(n/2)
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 = 1e06
# Plot the data
def Hantush_Jacob_fit(dlogs, dlogt, beta):
global T, S, B # 设置全局变量，值可以穿透程序
plt.style.use('default')
fig, ax = plt.subplots()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_aspect(1)
plt.xlabel('$\log 1/u$')
plt.ylabel('$W(u)$')
ax.grid(True)
T = Q/4/np.pi*np.float_power(10, dlogs) # 配线法不同，公式也不一样
S = 4*T/r**2*np.float_power(10, -dlogt) # 配线法不同，公式也不一样
u = r**2*S/4/T/x
u0 = r**2*S/4/T/t
ax.plot(x, hantush_jacob(u, beta),label="Hantush-Jacob曲线")
ax.plot(x, theis(u),label="Theis曲线")
ax.plot(t*np.float_power(10, dlogt), s*np.float_power(10, dlogs), '*', label="观测值")
plt.legend(
prop={'family': 'Simsun', 'size': 10}, handlelength=2,
title="图例",
title_fontproperties={'family': 'KaiTi', 'size': 12})
ax.set_title("配线法", fontproperties={'family': 'KaiTi', 'size': 12}) # 指定显示中文字体
plt.show()
RSS = np.sum((hantush_jacob(u0, beta) - s)**2)
print(' T(m^2/min): ', '{:.4f}'.format(T))
print(' S: ', '{:.4e}'.format(S))
if beta>0 :
print(' B: ', '{:.4e}'.format(r/beta))
print('Residual Sum of Squares: ', '{:.4e}'.format(RSS))
# ipywidgets 小部件控制参数, ipywidgets 有缓冲功能，同一个 Notebook 复制的代码得不到所期望的结果
dlogs_v=np.log10(T*4*np.pi/Q)
dlogt_v=np.log10(4*T/S/r/r)
widgets.interact(
Hantush_Jacob_fit,
dlogs = widgets.FloatSlider(
value=dlogs_v, min=dlogs_v-1.0, max=dlogs_v+1.0, step=0.05,
description=r"$\Delta\lg s$ [-]:", continuous_update=False,
readout_format='.2f', disabled=False),
dlogt = widgets.FloatSlider(
value=dlogt_v, min=dlogt_v-1.0, max=dlogt_v+1.0, step=0.05,
description=r'$\Delta\lg t$ [-]:', continuous_update=False,
readout_format='.2f', disabled=False),
beta = widgets.FloatSlider(
value=r/B, min=0.0, max=+3.0, step=0.01,
description=r'$\beta$ [-]:', continuous_update=False,
readout_format='.2f', disabled=False),
);
image.png
