直线图解法的程序实现（Python）

欢迎点评！如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。

需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions

用 pip 安装：

pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions

以下内容源自地下水动力学课程教学内容，程序可在 Jupyter notebook 中运行。

以 Jacob 公式的 配线为例

写成点斜式

通过调整 与 拟合直线，并用如下公式计算参数：

程序设计

程序设计流程为
库导入 数据准备 -> 绘图准备 widgets.interact 调整 与

直接上代码：

%matplotlib widget

import numpy as np
import math as math
import matplotlib.pyplot as plt
import ipywidgets as widgets
from wellfunction import *

# 控制小数的显示精度
np.set_printoptions(precision=4)

# 准备数据
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 = 0.0
ymax = math.ceil(max(s*10))/10
imin = math.floor(math.log10(min(t)))
imax = math.ceil(math.log10(max(t)))  
xmin = 10**imin
xmax =  10**imax
x=np.linspace(xmin, xmax, 50)

# 初始参数，最简单方法是取中间相邻两点用 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)

i0 = 0.183*Q/T
t0 = r**2*S/2.25/T

# Plot the data
def line_fit(a, b):  # a - 截距，b - 斜率
    global i0, t0, T, S    # 全局变量
    # 设置图形
    plt.style.use('default')
    fig, ax = plt.subplots(dpi=100)
    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.plot(t, s, '*', label="观测值")
    # 绘直线
    ax.plot(x, b*np.log10(x/a), label="拟合直线")
    # 计算调整后的参数
    T = 0.183*Q/b
    S = 2.25*T*a/r**2    
    
    plt.legend(
        prop={'family': 'Simsun'}, handlelength=6,
        loc=4,title="图例",
        title_fontproperties={'family': 'KaiTi'})
    # 指定图标题，显示中文
    ax.set_title("直线图解法", fontproperties={'family': 'KaiTi'})   

    plt.show()
    # 输出参数
    print('  T = {:.4f} m^2/min'.format(T))
    print('  S = {:.4e}'.format(S))

# 小部件互动
widgets.interact(
    line_fit,
    a = widgets.FloatSlider(
        value=t0, min=1, max=t0*5.0, step=.1,
        description='$t_0$ [-]:',
        continuous_update=False,
        readout_format='.1f',
        disabled=False),
    b = widgets.FloatSlider(
        value=i0, min=i0-0.1, max=i0+0.1, step=0.001,
        #value=i0, min=0.5*i0, max=2*i0, step=0.001,
        description='$slope$ [-]:',
        continuous_update=False,
        readout_format='.3f',
        disabled=False)
    );

image.png

尽管我们求出的 ，根据散点图的形态可以看出含水层可能有补给（参考：10. 降深 - 时间曲线与标准井函数的对比分析），因此还需要用其他方法继续分析水文条件。