Напишите функцию с NumPy для вычисления интеграла с определенным допуском - PullRequest
Новая
       85

Напишите функцию с NumPy для вычисления интеграла с определенным допуском

0 голосов
/ 26 апреля 2018

Я хочу написать пользовательскую функцию для числовой интеграции выражения (функции python или lambda) с определенным допуском.Я знаю, что с scipy.integrate.quad можно просто изменить epsabs, но я хочу написать саму функцию, используя numpy.

С этого поста Я знаю функцию: 

def integrate(f, a, b, N):
    x = np.linspace(a+(b-a)/(2*N), b-(b-a)/(2*N), N)
    fx = f(x)
    area = np.sum(fx)*(b-a)/N
    return area

дает мне численное интегрирование с N сегментами.Как я могу написать другую функцию или расширить ее, чтобы получить вход tol и увеличивать N, пока разница между двумя последующими вычислениями не станет меньше заданного допуска?

1 Ответ

0 голосов
/ 26 апреля 2018

Используя имеющуюся функцию, можно начинать с разумного N, например 5, и продолжать удваивать число, пока не будет достигнута требуемая точность.

def integrate_tol(f, a, b, tol):
    N = 5
    old_integral = integrate(f, a, b, N)
    while True:
        N *= 2
        new_integral = integrate(f, a, b, N)
        if np.abs(old_integral - new_integral) < tol:
            return (4*new_integral - old_integral)/3
        old_integral = new_integral

Простой тест: 

f = lambda x: np.exp(x)     
print(integrate_tol(f, -1, 1, 1e-9))
print(np.exp(1)-np.exp(-1))   # exact value for comparison

печать 

2.3504023872876028
2.3504023872876028

Нет никакой гарантии, что ошибка действительно меньше tol (но опять же, scipy.quad также не гарантирует этого).На практике ошибка будет намного меньше, чем tol, из-за использованного мною трюка, называемого Ричардсон экстраполяция : возвращаемое значение (4*new_integral - old_integral)/3 в целом намного точнее, чем сами новые или старые приближения.(Объяснение: поскольку integrate использует правило средней точки, каждое удвоение N уменьшает погрешность примерно в 4 раза. Таким образом, взятие комбинации 4*new_integral - old_integral почти сводит на нет остаточную ошибку в обоих этих результатах.)

(Примечание: в цикле while не рекомендуется начинать с N = 1; вероятно, этого будет недостаточно, и существует более высокий риск остановки слишком рано из-за некоторого числового совпадения, например, функция равна нулю в группемест.)

...