Оценка интегралов в Python 3.7;странное поведение - PullRequest
Новая
       23

Оценка интегралов в Python 3.7;странное поведение

0 голосов
/ 09 декабря 2018

Я написал некоторый код для аппроксимации интеграла:используя Python 3.7, но происходит странное поведение, которое дает мне неверный результат.Я вывел формулу следующим образом:

Пустьзатема такжедля n = 1,2,3 ...

Это реализовано в моем коде: 

import numpy as np

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

Это должно привести к I = 0,0000419426270488826, но мой код дает мне 50,40429353428721.Я пытался выяснить, что происходит, используя операторы print: 

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

Вы можете видеть, что вторая половина уравнения становится отрицательной на итерации 17, но я не понимаю, почему, так как обаЯ и п должны быть положительными.Я предполагаю, что именно здесь начинаются проблемы.Кто-нибудь знает, почему результат неверен, и если мое предположение верно?Я использую Mac OS X el Capitan 10.11.6

1 Ответ

0 голосов
/ 09 декабря 2018

Я думаю, это из-за ошибки переполнения.Если вы увеличите точность с плавающей точкой, используя Десятичное число , это решит проблему.Вот код: 

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))

Вот результат: 

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

Кстати, формулу необходимо исправить, заменив ln(x) на ln(2).

Обновлено: вероятно, название проблемы: потеря значимости

...