np.trapz of a function

Question

I'm having trouble figuring np.trapz out. I'm supposed to code a trapezoidal rule myself and then compare it to np.trapz However, there is a catch. Say the integral is from a to b. I'm supposed to find the integral for a=1 b=1, a=0 b=2, a=0 b=3... a=0 b=10 and plot these values. Here is my code for the trapezoidal function I made:

## trapezoidal ##
import numpy as np
import matplotlib.pyplot as plt

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))

data = [] ## data array

for b in range (0, 10):## for loop going through z2 as b

    ## definitions 
    a = 0
    N = 100
    h = (b-a)/N
    integral = 0

    integral = 0.5*(f(a) + f(b)) ## endpoints

    for i in range(1,N):
        integral += f(a + (i*h))
    integral = integral*h
    integral = integral/(1+b)
    data.append(integral) ## appending each iteration into data array

 print(data)

plt.plot([1,2,3,4,5,6,7,8,9,10],data)
plt.xlabel('z')
plt.ylabel('DA')

and here is what I've attempted for np.trapz though I think I'm hopelessly wrong. ## np.trapz ## import numpy as np import matplotlib.pyplot as plt

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))

x = np.arange(1,10,1)
##area = []

for i in range (0,10):
    x = [0,i]
    y = f/(1+i)
    area.append(np.trapz(y,x))

print(area)

plt.plot([1,2,3,4,5,6,7,8,9,10],area)
plt.xlabel('z')
plt.ylabel('DA')

Answer 1

I'll try to show how to implement the trapezoid rule, and compare it its implementation from numpy. So, you have a function def f(x) that you want to integrate from a to b.

You've not specified anything about error/tolerance acceptance but when comparing different methods for numerically evaluating integrals (.., equations or similar), it is often common to compare how their error scale when increased system size, number of iteration steps and so on.

The Trapezoid rule is a common text book method and there are tons of free literature (e.g. on Wikipedia). Hence, you easily verify your implementation's validity because its behavior is already well known. Also, it is always a good idea to test your method on a problem than can be solved analytically (by pen and paper). Doing so will quickly reveal errors in your code/implementation.

So let's do that! Consider the trivial function g(x) = x. Integrating g(x) from 0 to 10 is easily done and results in (1/2) * (10^2 - 0^2) = 50. Here's g(x) in python,

def g(x):
  return x

We now implement a function for evaluating the integral of any function func, using the trapezoid rule as defined on Wikipedia (link above),

def my_trapz(func, a, b, n_steps):
  X = np.linspace(a,b,n_steps)
  integral = (func(a)+func(b))/2.0 + sum([func(x) for x in X])
  return integral * (b-a)/n_steps

Here, func is a python function (def) that takes exactly one argument. a is the lower and b the upper bound of the integral, while n_steps is the number of x-values to evaluate within the range [a,b]. The larger n_steps is, the more accurate the integral will be. The numpy.linspace(a,b,n) function creates an array of n linearly spaced numbers between a and b. a we may now call,

a = 0.0
b = 10.0
n_steps = 100

integral = my_trapz(g, a, b, n_steps)
print(integral) # outputs 50.4999999..

You'll notice that the higher value you input for n_steps, the closer the result gets to 50, as expected. Now, we want to compare out implementation with the one in numpy. After reading its documentation, we see that the integral is calculated by,

X = np.linspace(a,b,n_steps)

integral = np.trapz([g(x) for x in X], dx=(b-a)/n_steps)
print(integral) # outputs 49.9023437.., slightly better than out implementation

Here, we've used a list comprehension directly into the function call. Alternatively, we could instead have done

g_values = []
for x in :
  g_values.append(g(x))
integral = np.trapz(g_values, dx=(b-a)/n_steps)

We can compare them, by agreeing on a set of different n_steps and investigate how the two methods perform: for each value of n_steps in the list n_steps_list, calculate the corresponding integral using both numpy and our own method. Then plot integral result as a function of n_steps. That is,

a = 0.0
b = 10.0
n_steps_list = [2**n for n in range(3,10)] # =[8, 16, 32, ..., 1024]

integral_np = []
integral_my = []

for n_steps in n_steps_list:
  X = np.linspace(a,b,n_steps)
  dx = (b-a)/n_steps
  integral_np.append(np.trapz([g(x) for x in X], dx=dx))
  integral_my.append(my_trapz(g, a, b, n_steps))

plt.scatter(n_steps_list, integral_np, color='g', label='np')
plt.scatter(n_steps_list, integral_my, color='r', label='my')
plt.xlabel('number of steps (resolution)')
plt.ylabel('Integral result')
plt.legend()
plt.show()

The resulting graph is shown below,

We see that both methods converge exponentially toward the expected value, 50. This is in agreement with the error-analysis provided by the Wikipedia article.

Now, I suggest you try to repeat this process, but by replacing my trivial g(x) with your original function f(x) and investigate what happens.

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))