使用Gauss-Laguerre求积法计算Gamma函数出错求助

你的代码结果偏差极大的核心原因是误用了高斯拉盖尔求积的被积函数形式，你的节点与权重匹配是正确的，主因出在被积函数的定义上。

错误根源 #

高斯拉盖尔求积的公式为：

∫₀^∞ e^(-x) f(x) dx ≈ Σ (w_i * f(x_i))

而Gamma函数的定义是：
Γ(n) = ∫₀^∞ x^(n-1) e^(-x) dx

对比可知，Gamma函数的积分正好对应高斯拉盖尔求积中f(x) = x^(n-1) 的场景，不需要在被积函数里额外乘以e^(-x)。你的代码里integrand(x)返回x**4 * math.exp(-x)，相当于计算的是∫₀^∞ x^4 e^(-2x) dx，完全偏离了Gamma函数的定义，导致结果错误。

修正后的代码 #

import math

# 正确的被积函数：仅保留x^(n-1)，无需乘e^(-x)
def integrand(x):
    return x**4

# 高斯拉盖尔求积的节点与权重（n=5，你的取值是正确的）
weights = [0.2635603197181409102030619, 1.413403059106516792218407, 3.59642577104072208112447, 7.08581000585883755692212, 12.64080084427578265943]
nodes = [0.1834346424956498, 0.9165752759150713, 2.4148302753646825, 4.2898664265427511, 7.209661773352202]

gamma_value = 0
for w, x in zip(weights, nodes):
    gamma_value += w * integrand(x)

print("Gamma(5)计算结果：", gamma_value)

验证结果 #

运行修正后的代码，输出结果会接近24（因数值积分精度，可能出现如23.999999999999996的微小误差），符合预期。

扩展说明 #

如果要计算任意n值的Gamma函数，只需修改integrand(x)中的指数为n-1，同时要确保使用的高斯拉盖尔节点和权重的阶数与求积点数一致。

内容的提问来源于stack exchange，提问作者Abishek S