Float overflow in scipy quadrature

Question

I am trying to calculate the following function:

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]

Answer 1

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

Answer 2

I just got some help from my professor, and he suggested to simplify the problem by assuming that

Which yields a simplified version of the integral

The two versions converge really quickly

I selected s_ == 100 as the point where I will use the simplified equation.

def style(r, gamma, upsilon, G):
    x_ = (r-R)/h
    B = upsilon/(2*G)
    def integrand(s_):
        if s_<100:
            return np.cos(s_*x_)/(((1+2*s_**2 + np.cosh(2*s_))/(np.sinh(2*s_) - 2*s_))*s_ + B/h*s_**2)
        else:
            return np.cos(s_*x_)/(s_ + B/h*s_**2)
    
                              
    int_out = quad_vec(integrand, 0, np.inf)
    return gamma/(2*np.pi*G) * int_out[0]

This works perfectly, just need to make sure to use float64, as float32 gets really jittery when the absolute value of x increases.