On the quad or trapz'd in ChemE heaven

| categories: | tags: | View Comments

What is the difference between quad and trapz? The short answer is that quad integrates functions (via a function handle) using numerical quadrature, and trapz performs integration of arrays of data using the trapezoid method.

Let us look at some examples. We consider the example of computing $$\int_0^2 x^3 dx$$. the analytical integral is $$1/4 x^4$$, so we know the integral evaluates to 16/4 = 4. This will be our benchmark for comparison to the numerical methods.

We use the scipy.integrate.quad command to evaluate this $$\int_0^2 x^3 dx$$.

from scipy.integrate import quad

ans, err = quad(lambda x: x**3, 0, 2)
print ans

4.0


you can also define a function for the integrand.

from scipy.integrate import quad

def integrand(x):
return x**3

ans, err = quad(integrand, 0, 2)
print ans

4.0


1 Numerical data integration

if we had numerical data like this, we use trapz to integrate it

import numpy as np

x = np.array([0, 0.5, 1, 1.5, 2])
y = x**3

i2 = np.trapz(y, x)

error = (i2 - 4)/4

print i2, error

4.25 0.0625


Note the integral of these vectors is greater than 4! You can see why here.

import numpy as np
import matplotlib.pyplot as plt
x = np.array([0, 0.5, 1, 1.5, 2])
y = x**3

x2 = np.linspace(0, 2)
y2 = x2**3

plt.plot(x, y, label='5 points')
plt.plot(x2, y2, label='50 points')
plt.legend()


The trapezoid method is overestimating the area significantly. With more points, we get much closer to the analytical value.

import numpy as np

x2 = np.linspace(0, 2, 100)
y2 = x2**3

print np.trapz(y2, x2)

4.00040812162


2 Combining numerical data with quad

You might want to combine numerical data with the quad function if you want to perform integrals easily. Let us say you are given this data:

x = [0 0.5 1 1.5 2]; y = [0 0.1250 1.0000 3.3750 8.0000];

and you want to integrate this from x = 0.25 to 1.75. We do not have data in those regions, so some interpolation is going to be needed. Here is one approach.

from scipy.interpolate import interp1d
import numpy as np

x = [0, 0.5, 1, 1.5, 2]
y = [0,    0.1250,    1.0000,    3.3750,    8.0000]

f = interp1d(x, y)

# numerical trapezoid method
xfine = np.linspace(0.25, 1.75)
yfine = f(xfine)
print np.trapz(yfine, xfine)

ans, err = quad(f, 0.25, 1.75)
print ans

2.53199187838
2.53125


These approaches are very similar, and both rely on linear interpolation. The second approach is simpler, and uses fewer lines of code.

3 Summary

trapz and quad are functions for getting integrals. Both can be used with numerical data if interpolation is used. The syntax for the quad and trapz function is different in scipy than in Matlab.

Finally, see this post for an example of solving an integral equation using quad and fsolve.

org-mode source

Integrating equations in python

| categories: | tags: | View Comments

A common need in engineering calculations is to integrate an equation over some range to determine the total change. For example, say we know the volumetric flow changes with time according to $$d\nu/dt = \alpha t$$, where $$\alpha = 1$$ L/min and we want to know how much liquid flows into a tank over 10 minutes if the volumetric flowrate is $$\nu_0 = 5$$ L/min at $$t=0$$. The answer to that question is the value of this integral: $$V = \int_0^{10} \nu_0 + \alpha t dt$$.

import scipy

nu0 = 5     # L/min
alpha = 1.0 # L/min
def integrand(t):
return nu0 + alpha * t

t0 = 0.0
tfinal = 10.0
V, estimated_error = quad(integrand, t0, tfinal)
print('{0:1.2f} L flowed into the tank over 10 minutes'.format(V))

100.00 L flowed into the tank over 10 minutes


That is all there is too it!

org-mode source

Integrating a batch reactor design equation

| categories: integration | tags: reaction engineering | View Comments

For a constant volume batch reactor where $$A \rightarrow B$$ at a rate of $$-r_A = k C_A^2$$, we derive the following design equation for the length of time required to achieve a particular level of conversion :

$$t(X) = \frac{1}{k C_{A0}} \int_{X=0}^X \frac{dX}{(1-X)^2}$$

if $$k = 10^{-3}$$ L/mol/s and $$C_{A0}$$ = 1 mol/L, estimate the time to achieve 90% conversion.

We could analytically solve the integral and evaluate it, but instead we will numerically evaluate it using scipy.integrate.quad. This function returns two values: the evaluated integral, and an estimate of the absolute error in the answer.

from scipy.integrate import quad

def integrand(X):
k = 1.0e-3
Ca0 = 1.0  # mol/L
return 1./(k*Ca0)*(1./(1-X)**2)

sol, abserr = quad(integrand, 0, 0.9)
print 't = {0} seconds ({1} hours)'.format(sol, sol/3600)
print 'Estimated absolute error = {0}'.format(abserr)

t = 9000.0 seconds (2.5 hours)
Estimated absolute error = 2.12203274482e-07


You can see the estimate error is very small compared to the solution.