Gauss-Legendre over intervals -x -> infinity：高效变换权重和节点的自适应算法

Question

好吧，我知道之前有人问过这个问题，并给出了一个缩放 [-1, 1] 间隔 [a, b] 的有限示例，但是没有人发布过如何将其概括为 [-a, Infinity] （正如下面所做的，但不是（还）快）。这也展示了如何使用多个实现调用复杂函数（无论如何在定量期权定价中）。有基准 quad 代码，后跟 leggauss，以及有关如何实现自适应算法的代码示例的链接。我已经解决了大多数链接的 adaptive algorithm 困难 - 它目前打印除法积分的总和以表明它可以正常工作。在这里您会找到将范围从 [-1, 1] 转换为 [0, 1] 到 [a, Infinity] 的函数（感谢@AlexisClarembeau）。要使用自适应算法，我必须创建另一个函数以将 [-1, 1] 转换为 [a, b]，然后反馈到 [a, Infinity] 函数。

import numpy as np
from scipy.stats import norm, lognorm
from scipy.integrate import quad

a = 0
degrees = 50

flag=-1.0000
F = 1.2075
K = 0.1251
vol = 0.43
T2 = 0.0411
T1 = 0.0047

def integrand(x, flag, F, K, vol, T2, T1):
    d1 = (np.log(x / (x+K)) + 0.5 * (vol**2) * (T2-T1)) / (vol * np.sqrt(T2 - T1))
    d2 = d1 - vol*np.sqrt(T2 - T1)
    mu = np.log(F) - 0.5 *vol **2 * T1
    sigma = vol * np.sqrt(T1)
    return lognorm.pdf(x, mu, sigma) * (flag * x*norm.cdf(flag * d1) - flag * (x+K)*norm.cdf(flag * d2))

def transform_integral_0_1_to_Infinity(x, a): 
    return integrand(a+(x/(1-x)), flag, F, K, vol, T2, T1) *(1/(1-x)**2); 

def transform_integral_negative1_1_to_0_1(x, a): 
    return 0.5 * transform_integral_0_1_to_Infinity((x+1)/2, a)

def transform_integral_negative1_1_to_a_b(x, w, a, b):
    return np.sum(w*(0.5 * transform_integral_0_1_to_Infinity(((x+1)/2*(b-a)+a), a)))

def adaptive_integration(x, w, a=-1, b=1, lastsplit=False, precision=1e-10):
    #split the integral in half assuming [-1, 1] range
    midpoint = (a+b)/2
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint)
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b)
    return interval1+interval2 #just shows this is correct for splitting the interval

def integrate(x, w, a): 
    return np.sum(w*transform_integral_negative1_1_to_0_1(x, a))

x, w = np.polynomial.legendre.leggauss(degrees) 
quadresult = quad(integrand, a, np.Inf, args=(flag, F, K, vol, T2, T1), epsabs=1e-1000)[0]
GL = integrate(x, w, a)
print("Adaptive Sum Result:")
print(adaptive_integration(x, w))
print("GL result"); 
print(GL)
print("QUAD result")
print(quadresult)

仍然需要用更少的维度提高速度和准确性，因为我无法手动调整 -a 的 degrees 范围以获得收敛。为了说明为什么这是一个问题，请改为输入这些值：a=-20、F=50，然后是运行。您可以增加 degrees=1000 并看到如果不智能地应用此 Gauss-Legendre 算法没有任何好处。我对速度的要求是每个循环达到 0.0004 秒，而我用 Cythonized 的最后一个算法大约需要 0.75 秒，这就是为什么我尝试使用高斯-勒让德的低度数、高精度算法。对于 Cython 和多线程，完全优化的 Python 实现的这一要求大约是每个循环 0.007s（非矢量化、循环缠身、低效的例程可能是每个循环 0.1s，degrees=20，即%timeit adaptive_integration(x,w).

我已经实施了一半的可能解决方案在第 5/6 页 http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals，adaptive integration 而间隔 a-b（在这种情况下，我写了 transform_integral_negative1_1_to_a_b function) where the interval divided in 2 (@0.5), the function then evaluated on these 1/2 intervals, and the sum of the two 0->0.5 + 0.5->1与整个范围 0->1 的函数结果进行比较。如果精度不在公差范围内，则范围进一步细分为 0.25 和 0.75，再次为每个子区间评估函数，并与之前的 1/2 区间总和 @0.5 进行比较.如果一侧在公差范围内（例如 abs(0->0.5 - (0->0.25 + 0.25->0.5)) < precision），但另一侧不在公差范围内，则拆分在公差范围内的一侧停止，但在另一侧继续，直到达到 precision。此时对区间的每一片的结果求和得到精度更高的全积分。

可能有更快更好的方法来解决这个问题。我不在乎，只要它又快又准。这是我遇到的集成例程的最佳描述，供参考http://orion.math.iastate.edu/keinert/computation_notes/chapter5.pdf奖励是100pts bounty + 15pts for answer acceptance。感谢您协助使此代码快速准确！

编辑：

这是我对 adaptive_integration 代码的更改 - 如果有人可以快速完成这项工作，我可以接受答案并奖励赏金。第 7 页 http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals 上的这个 Mathematica 代码执行我尝试过的例程。它在不能很好收敛的例程上工作，请参阅下面的变量。现在我的代码出错了：RecursionError: maximum recursion depth exceeded in comparison 在某些输入上，或者如果 degrees 设置得太高，或者当它工作时没有接近 quad 结果，所以有些东西这里显然是错误的。

def adaptive_integration(x, w, a, b, integralA2B, remainingIterations, firstIteration, precision=1e-9):
    #split the integral in half assuming [-1, 1] range
    if remainingIterations == 0:
        print('Adaptive integration failed on the interval',a,'->',b)
    if np.isnan(integralA2B): return np.nan

    midpoint = (a+b)/2
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint)
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b)
    if np.abs(integralA2B - (interval1 + interval2))  < precision : 
        return(interval1 + interval2)       
    else:
        return adaptive_integration(x, w, a, midpoint, interval1, (remainingIterations-1), False) + adaptive_integration(x, w, midpoint, b, interval2, (remainingIterations-1), False) 

#This example doesn't converge to Quad

# non-converging interval inputs
a = 0 # AND a = -250
degrees = 10

flag= 1
F = 50
K = 0.1251
vol = 0.43
T2 = 0.0411
T1 = 0.0047

print(adaptive_integration(x, w, -1, 1, GL, 500, False))

degrees=100的输出（计算GL和degrees=10000的输出以获得更好的初始估计，否则，算法显然总是与自己的精度一致并且不会调用每次都失败的自适应路径）：

GL result:
60.065205169286379
Adaptive Sum Result:
RecursionError: maximum recursion depth exceeded in comparison
QUAD result:
68.72069173210338

Answer 1

在 Titus A. Beu 的 "Numerical Programming: A Practical Guide for Scientists and Engineers Using Python and C/C++" 中，您可以在此处的代码示例 integral.py 和 specfunc.py 中找到方法：http://phys.ubbcluj.ro/~tbeu/INP/libraries.html 您调用函数 xGaussLag(a, deg) 从另一个 .py 文件调用 Laguerre 和 returns 你在 a 和 infinity 之间调整的 (x,w)。下面是如何设置它（注意上面 deg=80 它非常慢，我只是向你展示如何通过修改上面的行来应用它）：

x, w = np.array(xGaussLag(a,deg))
gauss = sum(w * integrand(x, flag, F, K, vol, T2, T1))

在 deg=80 上非常接近收敛（更快）但我只是将 eps=1e-13 放在 xGaussLag 中并用这些结果推送 deg=150，尽管如此比 quad 33%：

QUADPACK 解决方案：0.149221620346 错误：1.49870924498e-12 高斯勒让德解：0.149238273747 QUADPACK 和 Gauss-Legendre 的区别：1.66534003601e-05

在 Cython 中，这比直接 Python 快 6 倍顺便说一句，还是太慢了，所以我现在要尝试使用来自@Alexis 的答案的 "FastGL" 包，只是按照我的想法发布这对将来的其他 SO 用户很有用。

Answer 2

我认为代码可以完成工作：

import numpy as np 
import math

deg = 10
x, w = np.polynomial.legendre.leggauss(deg)

def function(x): 
    # the function to integrate
    return math.exp(-x)

def function2(x, a): 
    return function(a+x/(1-x))/((1-x)**2); 

def anotherOne(x, a): 
    return 0.5 * function2(x/2 + 1/2, a)

def integrate(deg, a): 
    sum = 0 
    x, w = np.polynomial.legendre.leggauss(deg)
    for i in range(deg): 
        print("sum({}) += {} * {} (eval in {})".format(sum, w[i], anotherOne(x[i], a), x[i]))
        sum += w[i]*anotherOne(x[i], a)
    return sum; 

print("result"); 
print(integrate(10, 1))

它结合了从 a 积分到 inf 的方程和改变积分边界的方程。

我希望它能解决你的问题（它至少对 exp(-x) 有效）:)

如果您想要内联计算，程序会执行以下总和： $enter image description here$

它是以下内容的组合：

并且：

$enter image description here$

并且：

Answer 3

跨越无限域的积分应该始终引起您的怀疑。毕竟，大多数 "simple" 功能甚至无法在那里集成！这就是为什么人们通常在其中使用 exp(-x) 或 exp(-x**2) 之类的阻尼项。果然，对于这两种情况，您有特定的集成规则：

Gauss-Laguerre 对于 exp(-x) 和
Gauss-Hermite 对于 exp(-x**2)。

因此，与其尝试将积分转换为有限域，不如尝试将它们转换为

这通常可以手工完成。

Gauss-Legendre over intervals -x -> infinity：高效变换权重和节点的自适应算法

Gauss-Legendre over intervals -x -> infinity: adaptive algorithm to transform weights and nodes efficiently

python

numpy

scipy

wolfram-mathematica

numerical-integration