使用梯形规则的指数积分函数的数值评估答案

【问题标题】：Numerical evaluation of the exponential integral function using trapezoidal rule使用梯形规则的指数积分函数的数值评估
【发布时间】：2020-06-02 01:44:47
【问题描述】：

我想用梯形法则对指数积分函数进行数值计算。这个函数定义为：

参考可用here。此函数已在某些库中可用，例如 scipy.special。由于某些原因，我不想使用这些库。相反，我需要通过梯形规则直接评估这个函数。我写了梯形规则并检查它以确保它工作正常。然后我用它来对 Ei 函数进行数值评估。不幸的是，结果不正确，例如我想评估等于 1.89511 的 Ei(1) 但我编写的代码返回 infinity 这是错误的。代码如下：

import numpy as np

# Integration using Trapezoidal rule
def trapezoidal(f, a, b, n):

    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

# Define integrand
def Ei(t):

    return - np.exp(-t) / t

# Define Ei(1)
A = trapezoidal(Ei, -1, np.inf, 20)

print (A)

# Checking Ei(1)
from scipy.special import expi
print (expi(1))

你知道我可以如何修改上面的代码以获得正确的结果吗？谢谢！

【问题讨论】：

标签： python-3.x math exponential integral numerical-integration

【解决方案1】：

1) 您不能定义以+inf 结尾的范围并将其分成 20 个部分。

相反，您可以选择任意权利限制并增加它，直到差异 abs(integral(limit(i+1))-integral(limit(i))) 变得可以忽略不计

2) 考虑零点的函数评估（如果发生）。它会导致除以零

如果参数太接近零，请尝试将其移动一点。

【讨论】：