On the quad, or trapz'd in ChemE heaven?

| categories: basic matlab | View Comments



On the quad, or trapz'd in ChemE heaven?

John Kitchin

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's 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.

clear all; clc; close all;

quad usage

This command is equivalent to $\int_0^2 x^3 dx$

i1 = quad(@(x) x.^3,0,2)
i1 =


we can also define the function handle separately.

f = @(x) x.^3
f(2) % should be 8
f = 


ans =


This command is equivalent to $\int_0^2 f(x) dx$

i1 = quad(f,0,2)
error = (i1 - 4)/4
% Note the error is very small!
i1 =


error =


Trapz usage

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

x = [0 0.5 1 1.5 2];
y = f(x);  % or you could say y = x.^3

i2 = trapz(x,y)
error = (i2 - 4)/4
i2 =


error =


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

figure; hold all
ezplot(f,[0 2])
legend 'numerical data' 'analytical function'
% The area of the trapezoids is greater than the area under the function
% due to the coarse sampling of $f(x)=x^3$.

with a finer sampling we have:

x = linspace(0,2); % 100 points between 0, 1
y = x.^3;
i3 = trapz(x,y)
i3 =


and a much smaller error.

error = (i3 - 4)/4
error =


Combining numerical data with quad

You might want to combine numerical data with the quad function if you want to perform integrals easily. Let's 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 don't have data in those regions, so some interpolation is going to be needed. Here is the analytical integral and by numerical quadrature

i_analytical = 1/4*(1.75^4-0.25^4)
i = quad(f,0.25,1.75)
% they are practically the same.
i_analytical =


i =


One way to do this with the data provided is to use interpolation on a fine-spaced x-grid, and then use trapz on the resulting interpolated vectors.

xfine = linspace(0.25,1.75);
yinterp = interp1(x,y,xfine);
i4 = trapz(xfine,yinterp)
error = (i4 - i_analytical)/i_analytical
i4 =


error =


Another approach is to create a function handle to do the interpolation for you. This combines numerical data with the quad function!

finterp = @(X) interp1(x,y,X);
i5 = quad(finterp,0.25,1.75)
error = (i5 - i_analytical)/i_analytical
i5 =


error =


these two approaches are actually identical, as they are both based on linear interpolation. The second approach is slightly more compact and easier to evaluate multiple integrals. Compare:

i6 = quad(finterp,0.28,1.64)
i6 =



xfine = linspace(0.28,1.64);
yinterp = interp1(x,y,xfine);
i7 = trapz(xfine,yinterp)
% One line of code vs. three.
i7 =



trapz and quad are functions for getting integrals. Both can be used with numerical data if interpolation is used.

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

% categories: basic matlab
% tags: integration
blog comments powered by Disqus