好的，我知道之前有人用一个有限的缩放示例问过这个问题 [-1, 1]间隔 [a, b] Different intervals for Gauss-Legendre quadrature in numpy但是没有人发布如何将其概括为 [-a, Infinity] (正如下面所做的，但不是(还)快)。这也展示了如何使用多个实现调用复杂函数(无论如何在定量期权定价中)。有基准quad代码，后跟 leggauss ，以及有关如何实现自适应算法的代码示例的链接。我已经完成了大部分链接 adaptive algorithm difficulties - 它目前打印除积分的总和以表明它工作正常。在这里您可以找到将范围从 [-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)

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

我已经实现了一半的可能解决方案在这里 http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals第 5/6 页，adaptive integration而区间 a-b (在这种情况下，我编写了 transform_integral_negative1_1_to_a_b 函数)其中间隔被分成 2 (@ 0.5 )，然后在这 1/2 间隔上评估函数，以及两个 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 )，但另一侧不在容差范围内，则 split 在容差范围内的一侧停止，但在另一侧继续，直到 precision到达了。此时将区间各切片的结果相加，得到精度更高的全积分。

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

编辑:

这是我对 adaptive_integration 的更改代码 - 如果有人可以快速完成这项工作，我可以接受答案并奖励赏金。第 7 页的 Mathematica 代码 http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals做我尝试过的例程。它在不能很好收敛的例程上工作，请参阅下面的变量。现在我的代码出错了: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

最佳答案

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

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) 有效):)

如果您想要内联计算，程序会执行以下总和:

它是以下的组合:

并且:

并且:

关于python - Gauss-Legendre 区间 -x -> 无穷大 : adaptive algorithm to transform weights and nodes efficiently，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/37367688/