在Python中求解双曲函数

Question

在水文学中，蒸发指数（EI）和干燥指数（DI）之间的关系用以下方程描述。

EI = (DI * tanh(1/DI) *(1+ sinh(DI) - scosh(DI))^1/2

但是，我们在 python 中很难解决这个问题。有没有一种我们不知道的简单方法可以做到这一点？

我们收到错误“无法在给定容差内找到根。(4.25389829408258985009 > 2.16840434497100886801e-19)”

我们找不到任何合适的猜测

将 sympy 导入为 sp

EI列表=[] DI 列表 = []

对于范围（6）内的 i：
EI = 平均值_ET[i]/平均值_perc[i]
EIlist.append(EI)
打印（EI）

DI = sp.symbols('DI') eq1 =sp.Eq(EI, (DI * sp.tanh(1/DI) *(1+ sp.sinh(DI) - sp.cosh(DI)))**1/2)
解决方案 = sp.nsolve(eq1, DI)

我们希望得到一个数值解

Answer 1

您可以使用标准牛顿-拉夫森方法通过反演来解决它：您的 EI 与 DI 的曲线似乎是单调递增的（如下图所示）。

鉴于平方根的出现，最容易从 EI^2 得到 DI。

import numpy as np
import matplotlib.pyplot as plt


def EIsq(D):                                                                      # Note: this is EI ** 2
    return D * np.tanh( 1.0 / D ) * ( 1.0 + np.sinh( D ) - np.cosh( D ) )


def EIsq_deriv( D ):
    return np.tanh( 1.0 / D ) * ( 1.0 + np.sinh( D ) - np.cosh( D ) )        \
         - 1.0 / D / np.cosh( 1.0 / D ) ** 2 * ( 1.0 + np.sinh( D ) - np.cosh( D ) )   \
         + D * np.tanh( 1.0 / D ) * (       np.cosh( D ) - np.sinh( D ) )


def NewtonRaphson( f, fprime, target, guess ):
    TOL = 1.0e-6
    x, xold = guess, guess + 1.0
    while abs( x - xold ) > TOL:
        xold = x
        x -= ( f( x ) - target ) / fprime( x )
    return x


# First plot the expected behaviour of EI against DI
DI = np.linspace( 0.001, 10.0, 1000 )
EI = np.sqrt( EIsq( DI ) )
plt.plot( DI, EI )

# Now choose some values of EI and try to work out which DI gave them
EImod = np.array( [ 0.1, 0.3, 0.5, 0.7, 0.8, 0.9, 0.95, 0.975, 0.985, 0.99 ] )
DImod = np.zeros( len( EImod ) )
for i in range( len( EImod ) ): DImod[i] = NewtonRaphson( EIsq, EIsq_deriv, EImod[i] ** 2, 0.5 )
plt.plot( DImod, EImod, 'o' )

plt.xlabel( "DI" )
plt.ylabel( "EI" )
plt.show()

输出：

Answer 2

这是一个有点有趣的令人讨厌的非线性方程，与双曲轨道不同，除了数值之外求解起来有些棘手。您需要有一个良好的初始猜测，否则数值求解器将会发散。

我采用 EI = y 和 DI = x，这样你的方程就变成：

y  = x.tanh(1/x)*sqrt(1 + sinh(x) - cosh(x))

this simplifies to

y  = x.tanh(1/x)*sqrt(1 + exp(-x))

很明显，对于您想要求解 x 的 y<0 which might make some solvers go AWOL. I'm assuming given a y >= 0，没有真正的解。

我的快速而肮脏的分析解决方案获得了一个可以使用牛顿拉夫森或更好的哈雷改进的启动器，如下所示。首先，我们忽略第二项，这是一个小的修正，因此求解更简单的近似值：

y ~= x.tanh(1/x)

quick sleight of hand using the Pade approx for tanh(1/x) gets us

y ~= 4x^2/(4x^2+1)

this can be inverted to give an answer within 50% of the real root

x ~= sqrt(y/(1-y))/2

这里最坏情况的错误大约是低两倍。可以通过乘以软糖系数 ~1.2x 来稍微改进。目的是有一些快速计算的东西，足以引导迭代数值方法，以便它能够收敛于根。如果这还不够好，那么我建议绘制该图并拟合一个有理多项式。

到目前为止，我快速浏览了一下，最佳效果如下：

y = (0.108209115 + 3.195883964x^2)/(2.31400472 + 3.159526242x^2)

解决 x 的问题应该会给你一个与正确答案相差百分之几的起始猜测。