CODEWITHSUNDEEP
python-3.xmathexponentialintegralnumerical-integration
Leo
Leo

Reputation: 509

Numerical evaluation of the exponential integral function using trapezoidal rule

I want to evaluate the exponential integral function numerically using trapezoidal rule.This function is defined as:

enter image description here

The reference is available here. This function is already available in some libraries for example scipy.special. For some reasons I do not want to use these libraries. Instead, I need to evaluate this function directly by the trapezoidal rule. I wrote the trapezoidal rule and checked it to make sure it works fine. Then I used it for numerical evaluation of the Ei function. Unfortunately the results is not correct for example I want to evaluate Ei(1) which is equal to 1.89511 but the code I have written returns infinity which is wrong. Here is the code:

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

Do you know how I can modify the above code so I can get correct results? Thanks!

Upvotes: 0

Views: 342

Answers (1)

MBo
MBo

Reputation: 80107

1) You cannot define range ending woth +inf and divide it onto 20 parts.

Instead you can choose arbitrary right limit and increase it until difference abs(integral(limit(i+1))-integral(limit(i))) becomes negligible

2) Consider function evaluation in zero point (if occurs). It causes division by zero

If argument is too close to zero, try to shift it a bit.

Upvotes: 1

Related Questions