pythonfloating-pointintegralcalculusquad

Float overflow in scipy quadrature


I am trying to calculate the following function:

formula

For every point x I need to evaluate the integral numerically since it can't be solved analytically. (At least sympy and WolframAlpha can't do it) I am using scipy.integrate.quad for this.
However, the sinh functions throw math range or float overflow errors when s is larger than ~300. I tried numpy functions and Python built-in math libraries.

Is there a way to properly calculate this?

Here is the equation implementation:

from scipy.integrate import quad
h = 0.0005
G = 1
R= 0.002
upsilon = 0.06
gamma = 0.072

def style_s(r, gamma, upsilon, G):
    x_ = (r-R)/h
    B = upsilon/(2*G) # substitute B for constants in last operand in denominator
    def integrand(s_):
        return np.cos(s_*x_)/(((1+2*s_**2 + np.cosh(2*s_))/(np.sinh(2*s_) - 2*s_))*s_ + B/h*s_**2)
    
    int_out = quad(integrand, 0, inf)
    return gamma/(2*np.pi*G) * int_out[0]

Solution

  • You can write

    sinh(2x) = 0.5(exp(2x)-exp(-2x)) = exp(2x)*0.5*(1-exp(-4x))

    cosh(2x) = 0.5(exp(2x)+exp(-2x)) = exp(2x)*0.5*(1+exp(-4x))

    Then multiply top and bottom of your fraction by exp(-2x) to get

    exp(-2x) sinh(2x) = 0.5*(1-exp(-4x))

    exp(-2x) cosh(2x) = 0.5*(1+exp(-4x))

    The decaying exponential doesn't seem to mind large s.

    Note that there's also an s/s-type behaviour near s=0, which you might want to remove. However, the integrator doesn't seem to worry about that.

    from math import pi
    import numpy as np
    from scipy.integrate import quad
    
    h = 0.0005
    G = 1
    R= 0.002
    upsilon = 0.06
    gamma = 0.072
    
    def style_s(r, gamma, upsilon, G):
        x = (r-R)/h
        B = upsilon/(2*G) 
    
        def integrand( s ):
            E2 = np.exp( -2 * s )
            E4 = np.exp( -4 * s )
            return h * (0.5*(1.0-E4) - 2*s*E2) * np.cos(s*x)             \
                       / s / ( h * 0.5*(1.0+E4) + B*s*0.5*(1.0-E4)  +   2*s**2*(h-B)*E2  + h*E2  )
        
        int_out = quad(integrand, 0, np.inf)
        return gamma/(2*pi*G) * int_out[0]
    
    r = 0.003
    print( style_s( r, gamma, upsilon, G ) )
    

    You don't give the value of r: I've guessed r=0.003

    I'm still getting a warning, but the output is

    8.023057532105406e-05