数值分析平时作业-牛顿迭代法最近 5 年炎热季节 10% 湿球温度及相应气象条件，是电厂冷却塔热力计算不可或缺的重要

这是我参与8月更文挑战的第12天，活动详情查看：8月更文挑战

最近 5 年炎热季节 10% 湿球温度及相应气象条件，是电厂冷却塔热力计算不可或缺的重要参数。我国许多气象观测站于 2002 年～ 2005 年间陆续停止了湿球温度的观测，给电厂冷却塔设计带来了一定的困难。

基于道尔顿蒸发定律和热量守恒原理，根据牛顿热传导公式等可以导出空气中的水汽压与干湿球温度、大气压力等的关系

用干、湿球温度求空气中水汽压的计算公式：

e = {E_t} - Ap({t_{\rm{d}}} - {t_{\rm{w}}})

{t_{\rm{w}}} = {t_{\rm{d}}} - \frac{{{E_t} - e}}{{Ap}}

{E_t} = e + Ap({t_{\rm{d}}} - {t_{\rm{w}}})

空气中相对湿度:

U = \frac{e}{{{E_{{t_{\rm{d}}}}}}} \times 100\%

式中td、tw分别为干、湿球温度(单位：℃)；e、p、Et分别为水汽压、气压、湿球温度对应的纯水平液(冰)面饱和水汽压(单位：hPa)；A为干湿表系数(单位：℃-1)，由干湿表类型、通风速度、湿球是否结冰而定；U为相对湿度(单位：%)；Etd为干球温度对应的纯水平液(冰)面饱和水汽压(单位：hPa)

纯水平液面饱和水汽压:

\begin{array}{l} {\rm{log}}({E_{\rm{t}}}) = 10.79574\left({1 - \frac{{{T_1}}}{T}} \right) - 5.028{\rm{log}}\left({\frac{T}{{{T_1}}}} \right) + \\ 1.50475 \times {10^{ - 4}}\left({1 - {{10}^{ - 8.2969\left({\frac{T}{{{T_1}}} - 1} \right)}}} \right) + 0.42873 \times \\ \;\;\;\;\;\;\;{10^{ - 3}} \times ({10^{4.70955\left({1 - \frac{T}{{{T_1}}}} \right)}} - 1) + 0.78614 \end{array}

由于 $Ap({t_{\rm{d}}} - {t_{\rm{w}}})$ 较小暂时忽略不计。可得等式：

{E_t} = U * \frac{{{E_{{t_{\rm{d}}}}}}}{100}

1、整体思路 ①输入干球温度、相对湿度

②求出干球温度求饱和蒸汽压（ Goff-Gratch方程直接算出）

③求出相对湿度求当前蒸汽压

④由当前蒸汽压倒推露点温度（牛顿迭代法倒推Goff-Gratch方程）

from sympy import diff, symbols, log


def f(x, y_0=0.0):
    """
    待求解的原函数
    :param x: float类型
    :param y_0: float类型。待求的x_0值，使得f(x_0)=y_0。
    :return: float类型
    """

    # Goff-Gratch方程，x=温度/K，y=饱和蒸汽压/hPa
    y = pow(10,
            (10.79574 * (1 - 273.16 / x)
             - 5.028 * log(x / 273.16, 10.0)
             + 1.50475 * 1e-4 * (1 - pow(10, -8.2969 * (x / 273.16 - 1)))
             + 0.4287 * 1e-3 * (pow(10, 4.76955 * (1 - 273.16 / x)) - 1)
             + 0.78614
             )
            )

    return y - y_0


def df(x):
    """
    原函数的导数
    :param x: float类型
    :return: float类型，返回df(x)。
    """

    t = symbols('t', real=True)  # 用于f求导的sympy.symbol类型
    y = diff(f(t), t, 1)  # 求f的导数，避免混淆以t做自变量

    return y.evalf(subs={t: x})  # 参数x代入自变量t求值


def newtonMethod(x_n, y_0, time=0):
    """
    牛顿迭代法
    :param x_n: float类型，本次迭代的x
    :param y_0: float类型，使f(x_n)值尽可能逼近y_0
    :param time: int类型，迭代次数显示
    :return: 迭代退出后返回x_n，使得f(x_n)=y_0
    """

    fx_n = f(x_n, y_0=y_0)  # 避免后续过程多次计算fx
    dfx_n = df(x_n)  # 避免后续过程多次计算dfx

    # print('第%s次迭代：\tx=%s,\tf(x) = %s,\tdf(x) =%s' %
    #       (str(time), str(x_n), str(fx_n), str(dfx_n),))

    if f(x_n, y_0=y_0) == 0.0:  # 如果f(x)恰好为0，则完成迭代
        # print('解为：\tx = %s,\tf(x) = %s' % (str(x_n), str(f(x_n))))
        return x_n
    else:
        x_n1 = x_n - fx_n / dfx_n

    dx_n1 = f(x_n1, y_0=y_0)
    if abs(fx_n - dx_n1) < 1e-6:  # 迭代完成条件
        # print('解为：\tx = %s,\tf(x) = %s' %
        #       (str(x_n1), str(dx_n1)))
        return x_n1  # 精度足够，返回解x_n1
    else:
        return newtonMethod(x_n1, y_0, time=time + 1)   # 精度未够，进入下次迭代


if __name__ == '__main__':
    dry_bulb = float(input('请输入干球温度/℃：'))
    relative_humidity = float(input('请输入相对湿度/%：'))
    print('________________')

    saturated_vapor_pressure = f(dry_bulb + 273.15)  # 饱和蒸汽压/hPa
    recent_vapor_pressure = saturated_vapor_pressure * relative_humidity / 100  # 当前蒸汽压/hPa
    absolute_humidity = 0.622 * relative_humidity / (1013.25 - relative_humidity)  # 绝对湿度  kg/kg绝干气。1013.25hPa是标准大气压
    print(recent_vapor_pressure)
    # 露点温度/℃，以牛顿迭代式求解（以干球温度为初始逼近值）
    dew_point_temperature = newtonMethod(dry_bulb + 273.15, recent_vapor_pressure) - 273.15

    print('绝对湿度：\t%s\t%s' % (round(absolute_humidity, 4), 'kg/kg绝干气'))
    print('露点温度：\t%s\t%s' % (round(dew_point_temperature, 1), '℃'))