log(x) 的泰勒级数

Question

我正在尝试评估自然对数 ln(x) 的泰勒多项式，在 Python 中以 a=1 为中心。我正在使用维基百科上给出的系列但是当我尝试像 ln(2.7) 这样的简单计算而不是给我接近 1 的东西时，它给了我一个巨大的数字。有什么明显的我做错了吗？

def log(x):
    n=1000
    s=0
    for i in range(1,n):
        s += ((-1)**(i+1))*((x-1)**i)/i
    return s

使用泰勒级数：

给出结果：

编辑：如果有人偶然发现，评估某个实数的自然对数的另一种方法是使用数值积分（例如，黎曼和、中点规则、梯形规则、辛普森规则等）来评估经常使用的积分定义自然对数；

Answer 1

该系列仅在 x 1，您将需要不同的系列。

比如这个（找到here）：

def ln(x): return 2*sum(((x-1)/(x+1))**i/i for i in range(1,100,2))

输出：

ln(2.7)        # 0.9932517730102833

math.log(2.7)  # 0.9932517730102834

请注意，随着 x 变大，收敛需要超过 100 个项（直到它变得不切实际）

您可以通过添加 x 的较小因子的对数来弥补这一点：

def ln(x):
    if x > 2: return ln(x/2) + ln(2)  # ln(x) = ln(x/2 * 2) = ln(x/2) + ln(2)
    return 2*sum(((x-1)/(x+1))**i/i for i in range(1,1000,2))

您也可以在基于 Taylor 的函数中执行此操作以支持 x>1：

def log(x):
    if x > 1: return log(x/2) - log(0.5) # ln(2) = -ln(1/2)
    n=1000
    s=0
    for i in range(1,n):
        s += ((-1)**(i+1))*((x-1)**i)/i
    return s

当 x 接近于零时，这些数列还需要更多项来收敛，因此您可能还希望在另一个方向上对它们进行处理，以保持计算的实际值在 0.5 和 1 之间：

def log(x):
    if x > 1:   return log(x/2) - log(0.5) # ln(x/2 * 2) = ln(x/2) + ln(2)
    if x < 0.5: return log(2*x) + log(0.5) # ln(x*2 / 2) = ln(x*2) - ln(2) 
    ...

如果性能是一个问题，您需要将 ln(2) 或 log(0.5) 存储在某处并重复使用，而不是在每次调用时都计算它

例如：

ln2 = None
def ln(x):
    if x <= 2:
        return 2*sum(((x-1)/(x+1))**i/i for i in range(1,10000,2))
    global ln2
    if ln2 is None: ln2 = ln(2)    
    n2 = 0
    while x>2: x,n2 = x/2,n2+1
    return ln2*n2 + ln(x)

Answer 2

程序是正确的，但Mercator series 有以下警告：

只要 -1

当x > 1 时，系列会发散，因此您不应期望结果接近 1。

Answer 3

python 函数math.frexp(x) 可用于修改问题，以便泰勒级数使用接近一的值。 math.frexp(x) 描述为：

将 x 的尾数和指数作为 (m, e) 对返回。 m 是一个浮点数 e 是一个整数，使得 x == m * 2**e 正好。如果 x 为零， 返回 (0.0, 0)，否则 0.5

使用math.frexp(x) 不应被视为“作弊”，因为它可能只是通过访问底层二进制浮点表示中的位字段来实现。不能绝对保证浮点数的表示将是 IEEE 754 binary64，但据我所知，每个平台都使用它。可以检查sys.float_info 以找出实际的表示细节。

很多like the other answer 可以使用如下标准对数恒等式吗：让m, e = math.frexp(x)。那么 log(x) = log(m * 2^e) = log(m) + e * log(2)。 log(2) 可以提前预计算到全精度，并且在程序中只是一个常数。这里有一些代码说明了计算 log(x) 的两个相似的泰勒级数近似值。每个系列中的术语数量是通过反复试验而不是严格分析确定的。

taylor1 实现 log(1 + x) = x¹ - (1/2) * x² + (1/3) * x^{3 ...}

taylor2 实现 log(x) = 2 * [t + (1/3) * t³ + (1/5) * t⁵ ...] , 其中 t = (x - 1) / (x + 1)。

import math
import struct

_LOG_OF_2 = 0.69314718055994530941723212145817656807550013436025

def taylor1(x):
    m, e = math.frexp(x)
    log_of_m = 0
    num_terms = 36
    sign = 1
    m_minus1_power = m - 1
    for k in range(1, num_terms + 1):
        log_of_m += sign * m_minus1_power / k
        sign = -sign
        m_minus1_power *= m - 1
    return log_of_m + e * _LOG_OF_2


def taylor2(x):
    m, e = math.frexp(x)
    num_terms = 12
    half_log_of_m = 0
    t = (m - 1) / (m + 1)
    t_squared = t * t
    t_power = t
    denominator = 1
    for k in range(num_terms):
        half_log_of_m += t_power / denominator
        denominator += 2
        t_power *= t_squared
    return 2 * half_log_of_m + e * _LOG_OF_2

这似乎适用于 log(x) 的大部分域，但是当 x 接近 1（并且 log(x) 接近 0）时，x = m * 2^e 提供的转换实际上会产生不太准确的结果。因此，更好的算法会首先检查 x 是否接近 1，例如 abs(x-1)