Численная оценка экспоненциальной интегральной функции с использованием правила трапеции - PullRequest
Новая
       81

Численная оценка экспоненциальной интегральной функции с использованием правила трапеции

0 голосов
/ 18 февраля 2020

Я хочу оценить экспоненциальную интегральную функцию численно, используя трапециевидное правило. Эта функция определяется как:

enter image description here

Ссылка доступна здесь . Эта функция уже доступна в некоторых библиотеках, например scipy.special . По некоторым причинам я не хочу использовать эти библиотеки. Вместо этого мне нужно оценить эту функцию непосредственно по правилу трапеции. Я написал правило трапеции и проверил его, чтобы убедиться, что оно работает нормально. Затем я использовал его для численной оценки функции Ei. К сожалению, результаты неверны, например, я хочу оценить Ei (1), равное 1,89511, но код, который я написал, возвращает бесконечность , что неверно. Вот код: 

import numpy as np

# Integration using Trapezoidal rule
def trapezoidal(f, a, b, n):

    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

# Define integrand
def Ei(t):

    return - np.exp(-t) / t

# Define Ei(1)
A = trapezoidal(Ei, -1, np.inf, 20)

print (A)

# Checking Ei(1)
from scipy.special import expi
print (expi(1))

Знаете ли вы, как я могу изменить вышеуказанный код, чтобы получить правильные результаты? Спасибо!

1 Ответ

1 голос
/ 18 февраля 2020

1) Вы не можете определить окончание диапазона +inf и разделить его на 20 частей.

Вместо этого вы можете выбрать произвольный правильный предел и увеличивать его до тех пор, пока разница abs(integral(limit(i+1))-integral(limit(i))) не станет незначительной

2) Рассмотреть оценку функции в нулевой точке (если происходит). Это вызывает деление на ноль

Если аргумент слишком близок к нулю, попробуйте сместить его немного.

...