Monday, June 13, 2011

Calculating the Area Under a Surface

starLoren on the Art of MATLAB
June 13, 2011 9:25 AM
by Loren

Calculating the Area Under a Surface

Recently there was an email making the rounds at MathWorks about how to calculate the area under a surface. Not surprisingly, there were several methods chosen, based on each sender's proclivities. Here are some of the ways.

Contents

Data Already on a Regular Grid

If you have data already on a regular grid, you can simply call trapz twice, once along the X dimension, and once along Y. Here's an example. Let's first make some random x and y points.

xdata = [0; rand(100,1); 1]; ydata = [0; rand(100,1); 1]; x = sort(xdata); y = sort(ydata);

For now, let's calculate our gridded sampled function. Then use trapz twice.

[X,Y] = meshgrid(x,y); Z = X.^2.*sin(3*(X-Y)); trapz(y,trapz(x,Z,2),1)
ans =       0.13173 

Let's compare the result to one using dblquad.

F = @(x,y)(x.^2).*sin(3*(x-y)); dblquad(F,0,1,0,1)
ans =       0.13173 

Scattered Data : Finding the Convex Hull

MATLAB has the ability to deal with scattered data in a variety of ways. I've just shown, using meshgrid a technique for integrating a function that you know by sampling, and comparing this result to numerically integrating the same function that generated the samples. I'll take advantage of some of the newer computational geometry functionality in MATLAB in this next attempt.

load seamount minz = min(z); zadj = z-min(z);

Note that I "normalized" the depth data (z), since it is negative, and changed it into the array zadj that is all non-negative. Next I create an interpolant from the scattered data. Then integrate this function using quad2d, but not before trimming the minimum and maximum values for x and y so I won't run into extrapolated values (which are NaN).

F2 = TriScatteredInterp(x,y,zadj); q1 = quad2d(@(x,y) F2(x,y),211.1,211.4,-48.35,-48,'AbsTol',0.01)
q1 =        136.04 

Scattered Data : Finding the Area

If you have access to Curve Fitting Toolbox, you can take advantage of the relatively new capability for fitting surfaces. I'll use the same seamount data as before, created a fit to the data, and, from the function that comes from the fit, use quad2d again to compute the area.

f = fit([x, y], zadj, 'linearinterp'); q2 = quad2d(f,211.1,211.4,-48.35,-48,'AbsTol',0.01)
q2 =        136.04 

How Do You Work with Scattered Data?

If you have scattered data representing a surface, how do you think about it and work with it - numerically or by fitting a function of some sort? Do you have other techniques? Let me know here.


Get the MATLAB code (requires JavaScript)

Published with MATLAB® 7.12

Computational Geometry How To Interpolation & Fitting


Dr. Art Trembanis
CSHEL
109 Penny Hall
Department of Geological Sciences
The College of Earth, Ocean, and Environment
University of Delaware
Newark DE 19716
302-831-2498

"Education is not the filling of a pot, but the lighting of a fire." -W. B. Yeats

No comments: