MATH3315 -Project 2_Ocean Volume Estimation - Solved

Numerical integration of a function of one variable is called quadrature. The analogous situation of integrating a function of two variables is sometimes called cubature. In this project you will explore the extension of the trapezoid formula to this situation, and apply it to find the volume of water in a patch of the ocean.

Consider the problem RRR f (x, y)dx d y , where R is the rectangle [x0,xm]×[y0, yn]. We discretize each variable using the space steps h and k, respectively, so that xi = x0 + ih for i = 0,...,n, and yj = y0 + jk for j = 0,...,m. (This requires that h = (xm − x0)/m and k = (yn − y0)/n.) The result is a decomposition of R into rectangles with the nodes at the corners. The strategy is to compute a volume over each little rectangle and sum them over the whole domain. Here’s an example of how the discretization process works in MATLAB:

Construct the domain of integration.

xLeft = 0; xRight = 4; yLeft = -1; yRight = 1; fill([xLeft xRight xRight xLeft],[yLeft yLeft yRight yRight],0.8*[1 1 1]) axis equal, hold on

 

These are the individual variable discretizations.

m                   = 6; h = (xRight-xLeft)/m;x = xLeft + h*(0:m)’;

n                     = 5; eta = (yRight-yLeft)/n;y = yLeft + eta*(0:n)’;

set(gca,’xtick’,x,’ytick’,y)

2

 

Now we construct a 2D grid out of them.

[X,Y] = ndgrid(x,y); plot(X,Y,’b.’) plot(X,Y,’k--’)

 

                      0.6667    1.3333        2        2.6667    3.3333        4

Both X and Y are (m +1)×(n +1) matrices, and satisfy the identity (xi, yj)=(Xi j,Yi j).

point = [x(3),y(2)] pointAgain = [ X(3,2),Y(3,2) ]

point =

1.3333
-0.6000
pointAgain =

1.3333
-0.6000
At the corners of each rectangle we have values of z according to Zi j = f (xi, yj). We can model the surface z = f (x, y) over one of these rectangles using the function Li j(x, y)=a +b(x −xi)+c(y −yj)+ d(x −xi)(y −yj), where a,b,c,d are determined by interpolation conditions at the four corners of the rectangle. (These coefficients depend on i and j, but that is left out to make the notation clearer.) We get a solid defined over the rectangle looking like this:

Here is a single rectangle.

h = 0.1; x = [0 h]’; eta = 0.2; y = [-1 -1+eta]’;

[X,Y] = ndgrid(x,y)

X =

0
0
0.1000
0.1000
Y =

-1.0000
-0.8000
-1.0000
-0.8000
These are made-up values at the corners of the rectangle.

Z = [ 1.1, 0.9; 0.75, 2 ];

Here are coefficients of the interpolant.

A = [ ones(4,1), X(:), Y(:), X(:).*Y(:) ]; abcd = A \ Z(:)

abcd =

0.1000 69.0000 -1.0000

72.5000
And here is the resulting surface and the outlines of the solid. (Don’t worry about the fancy plotting commands here.)

ezsurf(@(x,y) abcd(1) + abcd(2)*x + abcd(3)*y + abcd(4)*x.*y, [x;y] ) hold on, shading interp plot3(X,Y,0*Z,’k’,’linew’,2) plot3([1;1]*X(:)’,[1;1]*Y(:)’,[zeros(1,4);Z(:)’],’k’,’linew’,2) axis tight



Computational content:

1.    Download and install the GeoMapApp at http://www.geomapapp.org.

2.    Starttheappandagreetothelicense. Youshouldgetawindowwithaprojectedmapoftheworld, in latitude-longitude corrdinates as measured in degrees of arc, and some GUI widgets. Use the magnifying glass tool to select a small rectangle over the ocean. You should try to get it to include about 1-2 degrees of latitiude (the y axis).

3.    Click on the icon that looks like a grid. After a few seconds another window should open with a histogram in it. Make sure the dropdown menu at the top of this window says “GMRT Grid” with some version number. Click on the floppy disk icon, select “Grid: .XYZ (ascii format)”, and save the file to your computer. At this point you can close GeoMapApp.

4.    You can import this data into MATLAB by double-clicking on it in the “Current Folder” tab of the MATLAB Desktop, and using the resulting dialog. What you will get is three columns of x, y , and z data on a regular rectangular grid, where x and y are reported in degrees of arc, and z will be negative to represent depths below sea level.

5.    Afterconvertingallthedataintokmunits, use it and the formula you derived above to estimate the amount of water in your patch of ocean.