CODEWITHSUNDEEP
pythonpython-3.xmathcalculus
tail_recursion
tail_recursion

Reputation: 740

Taylor series for log(x)

I'm trying to evaluate a Taylor polynomial for the natural logarithm, ln(x), centred at a=1 in Python. I'm using the series given on Wikipedia however when I try a simple calculation like ln(2.7) instead of giving me something close to 1 it gives me a gigantic number. Is there something obvious that I'm doing wrong?

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

Using the Taylor series:

enter image description here

Gives the result: enter image description here

EDIT: If anyone stumbles across this an alternative way to evaluate the natural logarithm of some real number is to use numerical integration (e.g. Riemann sum, midpoint rule, trapezoid rule, Simpson's rule etc) to evaluate the integral that is often used to define the natural logarithm;

enter image description here

Upvotes: 2

Views: 4830

Answers (4)

JoEmmyA
JoEmmyA

Reputation: 1

My answer is just using the Taylor series for In(x). I really hope this helps. It is simple and straight to the point. enter image description here

Upvotes: -2

President James K. Polk
President James K. Polk

Reputation: 42019

The python function math.frexp(x) can be used to advantage here to modify the problem so that the taylor series is working with a value close to one. math.frexp(x) is described as:

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.

Using math.frexp(x) should not be regarded as "cheating" because it is presumably implemented just by accessing the bit fields in the underlying binary floating point representation. It isn't absolutely guaranteed that the representation of floats will be IEEE 754 binary64, but as far as I know every platform uses this. sys.float_info can be examined to find out the actual representation details.

Much like the other answer does you can use the standard logarithmic identities as follows: Let m, e = math.frexp(x). Then log(x) = log(m * 2e) = log(m) + e * log(2). log(2) can be precomputed to full precision ahead of time and is just a constant in the program. Here is some code illustrating this to compute the two similar taylor series approximations to log(x). The number of terms in each series was determined by trial and error rather than rigorous analysis.

taylor1 implements log(1 + x) = x1 - (1/2) * x2 + (1/3) * x3 ...

taylor2 implements log(x) = 2 * [t + (1/3) * t3 + (1/5) * t5 ...], where 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

This seems to work well over most of the domain of log(x), but as x approaches 1 (and log(x) approaches 0) the transformation provided by x = m * 2e actually produces a less accurate result. So a better algorithm would first check if x is close to 1, say abs(x-1) < .5, and if so the just compute the taylor series approximation directly on x.

Upvotes: -1

Alain T.
Alain T.

Reputation: 42129

That series is only valid when x is <= 1. For x>1 you will need a different series.

For example this one (found here):

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

output:

ln(2.7)        # 0.9932517730102833

math.log(2.7)  # 0.9932517730102834

Note that it takes a lot more than 100 terms to converge as x gets bigger (up to a point where it'll become impractical)

You can compensate for that by adding the logarithms of smaller factors of 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))

which is something you can also do in your Taylor based function to support 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

These series also take more terms to converge when x gets closer to zero so you may want to work them in the other direction as well to keep the actual value to compute between 0.5 and 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) 
    ...

If performance is an issue, you'll want to store ln(2) or log(0.5) somewhere and reuse it instead of computing it on every call

for example:

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)

Upvotes: 7

PaSTE
PaSTE

Reputation: 4633

The program is correct, but the Mercator series has the following caveat:

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

The series diverges when x > 1, so you shouldn't expect a result close to 1.

Upvotes: 3

Related Questions