CODEWITHSUNDEEP
python-3.xfor-loopmathintegral
B.Ing
B.Ing

Reputation: 67

Evaluating integrals in Python 3.7; strange behaviour

I written some code to approximate an integral: using Python 3.7, but some strange behaviour is happening which is giving me the wrong result. I derived the formula in the following way:

Let then and for n = 1,2,3...

This implemented in my code:

import numpy as np

I = 1
for n in range(1,21):
    I = 2*(np.log(2))**n - n*I

This should result in I = 0.0000419426270488826, but my code is giving me 50.40429353428721. I've been trying to figure out what is happing using print statements:

print("Iteration: ",n)
print("first half: ",2*(np.log(2))**n)
print("second half: ", n*I)
print("New I: ",I, "\n")

You can see that the second half of the equation turns negative in iteration 17, but I don't understand why, since both I and n should be positive. My guess is that this is where the problems start to occur. Does anyone know why the result is incorrect, and if my assumption is correct? I'm using Mac OS X el Capitan 10.11.6

Upvotes: 0

Views: 112

Answers (1)

lenhhoxung
lenhhoxung

Reputation: 2766

I think it's because of the overflow error. If you increase the float accuracy by using Decimal, it will solve the problem. Here is the code:

import numpy as np
from decimal import *
getcontext()

n = 20
print('n=',n)
def In(n):
    if n==0:
        return 1
    else:
        return 2*Decimal(2).ln()**n - n*In(n-1)

print('I{}={}'.format(n,float(In(n))))

# check with an existing function of Scipy
import scipy.integrate as integrate
result = integrate.quad(lambda x: np.log(x)**n, 1, 2)
print('I{}={}'.format(n,result))

Here is the output:

n= 20
I20=4.1942535120159466e-05
I20=(4.1942627048882645e-05, 7.29235597161084e-18)

By the way, the formula needs to be corrected by replacing ln(x) with ln(2).

Updated: Probably the name of the problem is loss of significance

Upvotes: 3

Related Questions