二重指数型積分公式とは

Double exponential formula (DE公式), tanh-sinh quadratureとも。
Wikipediaや提案者森さんのPDF参照。数値積分をする上で、端点での特異性が強い関数の場合に強い（し、普通の関数に対して使ってもいい感じ）。
scipy.integrate.quadだと端点で特異性が強いと端点での特異性に引きずられて計算時間がかかるわ精度も当てにならないわで大変な場合があるが、
この公式を使って変数変換してしまえばあとは台形公式を使うだけで計算量もそこそこにいい感じの精度が出る（というか台形公式を使ったときに最も効率よく精度が出るような仕掛けになっている）。
なんかMathematicaでも実装されているらしいけどMathematicaは高いし貧乏人的には使いたくないので、pythonで実装してみる。

ライブラリ・モジュール

Python: quadpy
C++: Boost
Mathematica: NIntegrateでMethod->"DoubleExponential"と指定

quadpyがあるならそれ使えばいいという話もあるが、今回は似たような積分をパラメータだけ変えて大量に積分する必要があったので、numpyのbroadcastを使って全部いっぺんに求めたい（単純にForで回すとPythonなので遅い可能性を危惧して）。

quadpy is fully vectorized, so if you like to compute the integral of a function on many domains at once, you can provide them all in one integrate() call, e.g.,...

ということなので、実装する必要はなかった。ちゃんちゃん。

実装

~~半無限区間だけ必要だったのでとりあえず雑な実装をしてみた：~~
2019/10/16: 重み計算用の関数を分離、簡単なDocstringを追加

dequad.py

from numpy import sinh,cosh,exp,log,pi,arange,isnan,isinf,float64,float128
import numpy as np
import matplotlib.pyplot as plt

DEBUG = False
BUF_DEBUG = None

@lru_cache(maxsize = None)
def generate_x_w(a,b,width,mN,pN,xp=np):
    '''
    x = phi(t)
    w = phi'(t)
    '''
    #print("generate_x_w in",a,b,width,mN,pN)
    pi2 = xp.pi/2
    ts = width*xp.arange(-mN,pN)#.astype(dtype)
    if np.all(not np.isinf(a)) and np.all(np.isinf(b)):
        ps = pi2*xp.sinh(ts)
        xs = xp.exp(ps) + a
        ws = width *pi2 *xp.cosh(ts)*xp.exp(ps)    
        
    elif np.all(np.isinf(a)) and np.all(np.isinf(b)):
        ps = pi2*xp.sinh(ts)
        xs = xp.sinh(ps)
        ws = width *pi2 *xp.cosh(ts)*xp.cosh(ps)  
        
    elif np.all(not np.isinf(a)) and np.all(not np.isinf(b)):
        ps = pi2*xp.sinh(ts)
        xs = (b-a)/2*xp.tanh(ps) + (a+b)/2
        ws = width * (b-a)/2 *pi2 *xp.cosh(ts)/(xp.cosh(ps)**2)
        
    return xs,ws

def dequad(func,a,b,width=5e-3,pN=1000,mN=1000,axis=None,dtype=float64,show_fig=False,show_integrand_array=False,ignore_nan=False,xp=np):
    '''
    Double exponential quadrature.

    func: func(ndarray_in) = ndarray_out
    axis: define the axis of ndarray_out to use in the integration.
        Example: 
            For functions like 
                func: array(n,) -> array(m,n),
            we can obtain the integral of "func" along the the last axis (1,n):
                f(x) * w ~ (m,n) * (n,) ~ broadcast~ (m,n) * (1,n)
            by "axis=1". 
    '''
    
    #print(a,b,width,mN,pN)
    #print(generate_x_w(a,b,width,mN,pN,xp=np))
    xs,ws = generate_x_w(a,b,width,mN,pN,xp=np)
    wsfs = ws*func(xs)
    isnan = xp.isnan(wsfs)
    isinf = xp.isinf(wsfs)
    wsfs[isnan | isinf] = 0
    
    print(wsfs[wsfs!=0]) if DEBUG else None
    
    if np.any(isnan) and not ignore_nan:
        raise TypeError("dequad_hinf: wsfs is nan! {}".format(wsfs))
    if show_fig:
        if len(wsfs.shape)>1:
            np.set_printoptions(threshold=20)
            display(wsfs) if show_integrand_array else None
            plt.plot(width*xp.arange(-mN,pN),wsfs.T)
            plt.legend(wsfs.sum(axis=axis))
        else:
            plt.plot(ts,wsfs)
    return (wsfs).sum(axis=axis)

基本的にnumpy.arrayを使うようにと、端点での指数の発散をどうにかできればよし。

時間のあるときにもうちょっとまともに修正したい。