Trying to Understand Numerical Integration Using Python Syntax

Question

I am trying to learn to program in Python. While I will use the built-in functions in the future, I want to understand how to program the routines from scratch as I believe it will better help me understand Python syntax and coding much better. I have been trying to write a program to numerically integrate a function by direct summation. I have been unable to figure out what changes I need to make to my program to get this to work in Python and some help would be appreciated.

from math import *
import numpy as np

N = int(input("Enter value for how many sums you wish to perform: "))
a = float(input("Enter the lower integration bound: "))
b = float(input("Enter the higher integration bound: "))

def Integrate(N, a, b):
    def f(x):
        return x**2 * e**(-x)*np.sin(x)
    val = 0
    out = 0
    for i in range(N):
        function += f(a+((i-(1/2))*((b-a)/N)))
    out = ((b-a)/N)*val
    return out

I tried the above code and I am not even obtaining an output, so I suspect there is something regarding the Python language I am missing.

Answer 1

Basic example of numerical integration, Riemann integral with a uniformly partitioned interval.

For convenience the function is passed as argument.

def Integrate(N, a, b, f):
    # Riemann integral, evaluation at the midpoint of each interval
    width = (b-a) / float(N)
    return sum(f(a+(i+.5)*width)*width for i in range(N))


N = 10
a = 0
b = 1

# testing the area of a triangle

# test function - diagonal
def diagonal(x): return x

numerical_approx = Integrate(N, a, b, f=diagonal)
exact_area = (b-a)*diagonal(b) / 2 # base * height / 2

error = abs(numerical_approx - exact_area
print(f"{error):.20f}")
#0.00000000000000011102

The fact that the difference of areas is not zero is due to computational limitations of the float data types.

---

For the original problem

def f(x):
    return x**2 * np.e**(-x)*np.sin(x)

N, a, b = 10, 0, 1

I = Integrate(N=N, a=a, b=b, f=f)
print(I)
#0.10225500155762644

for a rough check you can use numpy.trap which use a trapezoidal approximation rule. Also if the result could be slightly different it is a good indicator (any integrable function will converge to the same value independently from the partion)

width = (b-a) / float(N)
y = [f(a+(i)*width) for i in range(N+1) ]
I_trap = np.trapz(y, dx=width)
print(I_trap)
#0.10289249898529564