Solving hyperbolic function in python

Question

In hydrology, the relationship between evaporative index (EI) and dryness index (DI) is described with the following equation.

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

However, we have trouble solving this in python. Is there an easy way to do this we're not aware about?

We get the error "Could not find root within given tolerance. (4.25389829408258985009 > 2.16840434497100886801e-19)"

we cant find any guesses which fits

import sympy as sp

EIlist = []
DIlist = []

for i in range(6):   
  EI = mean_ET[i]/mean_perc[i]   
  EIlist.append(EI)   
  print(EI)

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

We are hoping to get a numerical solution.

Answer 1

You can solve it by inversion using the standard Newton-Raphson method: your curve of EI against DI appears to be monotonically increasing (as the graph below shows).

Given the appearance of the square root it is easiest to get DI from EI^2.

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()

Output:

Answer 2

It is a faintly interesting nasty non-linear equation not unlike hyperbolic orbits and somewhat tricky to solve other than numerically. You need to have a good initial starting guess or a numerical solver will diverge.

I'm taking EI = y and DI = x so that your equation becomes:

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

this simplifies to

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

It is immediately apparent that there are no real solutions for y<0 which might make some solvers go AWOL. I'm assuming given a y >= 0 you want to solve for x.

My quick and dirty analytical solution to get a starter that you can refine using Newton-Raphson or better Halley is as follows. First we ignore that second term which is a small correction so solving the simpler approx:

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

Worst case error here is about a factor of two low. Could be improved a bit by multiplying through by a fudge factor ~1.2x. The aim is to have something quick to compute that is good enough to bootstrap the iterative numerical method so that it will be able to converge on the root. If that isn't quite good enough then I suggest plot the thing and fit a rational polynomial to it.

Having come this far I had a quick look and the optimum looks something like:

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

Solving that for x should give you a starting guess within a few percent of the right answer.

Answer 3

Note, first, that you have "x**1/2" (where x is part of your expression) and this just means x/2. Either use sqrt or write x**(1/2).

Equations that are ill-behaved in one form can be better behaved in another. First simplify the function:

>>> from sympy import *
>>> EI,DI = var("EI DI")
>>> eq = EI - sqrt(DI*(sinh(DI) - cosh(DI) + 1)*tanh(1/DI))
>>> seq = simplify(eq); seq
EI - sqrt(DI*(1 - exp(-DI))*tanh(1/DI))

Then solve for one of the 3 DI by replacing the other forms with dummy variables while doing so as described here. Two alternate forms are then

e = [
Eq(DI, EI**2*exp(DI)/((exp(DI) - 1)*tanh(1/DI))),
Eq(DI, log(DI*tanh(1/DI)/(DI*tanh(1/DI) - EI**2)))]

The first form behaves well for a wide range of evaporative indices

>>> [nsolve(e[0].subs(EI,j),1) for j in (.001,.1,.2,.5,.7,.99)]
[0.00100, 0.103, 0.211, 0.597, 1.02, 4.99]