log(x) 的泰勒级数

Question

我正在尝试计算以 Python 中的 a=1 为中心的自然对数 ln(x) 的泰勒多项式。我正在使用维基百科上给出的系列，但是当我尝试像 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

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

The series converges to the natural logarithm (shifted by 1) whenever −1 < x ≤ 1.

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

Answer 2

该系列仅在 x <= 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))

您也可以在基于泰勒的函数中执行此操作以支持 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 3

python 函数 math.frexp(x) 在这里可以用来修改问题，使泰勒级数的值接近 1。 math.frexp(x) 描述为：

Return the mantissa and exponent of x as the pair (m, e). m is a float and e is an integer such that x == m * 2**e exactly. If x is zero, returns (0.0, 0), otherwise 0.5 <= abs(m) < 1. This is used to “pick apart” the internal representation of a float in a portable way.

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

很多你可以使用标准对数恒等式如下：设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³ ...

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) < .5，如果是，则直接在 x 上计算泰勒级数近似值。

log(x) 的泰勒级数

Taylor series for log(x)

python

math

calculus

python-3.x