6.54 area problem (9.11.98)

6.54.1 Powerhouse Museum

How do I calculate the area under a curve in Maple when I haven’t got a formula but experimental data?

Please send response to newsgroup or directly to myself: [e-mail] (omit the nospam for the true address).

6.54.2 John Little (11.11.98)

There are any number of standard techniques from numerical analysis for this kind of calculation – they are called numerical quadrature rules. A simple, commonly used example of these methods is Simpson’s Rule. If you have an odd number n = 2k+1 of data points

          (x_1,f(x_1)), (x_2,f(x_2)), ... , (x_n, f(x_n)) 
with x_0 < x_1 < ... < x_n and with x_i - x_{i-1} = h for all i,
 

(i.e. the independent variable values are *equally spaced*) then the Simpson’s Rule approximation to the integral of f over the interval \([x_1,x_n]\) (the "area under the curve" that you want) is

  (f(x_1) + 4f(x_2) + 2f(x_3) + 4f(x_4) + ... + 4f(x_{n-1}) + f(x_n))*h/3
 

(the 4’s and 2’s alternate in the summation). If you have the dependent variable values stored in a 1-dimensional array (or list) you can compute this in Maple by the following procedure (assuming \(n = 2k + 1\) is odd):

  ListSimpson := proc(fvals,n,h) 
 
     local Simpson,j; 
 
     Simpson := fvals[1] + fvals[n]; 
     for j from 2 to n-1 do 
       if type(j,even) 
         Simpson := Simpson + 4*fvals[j] 
        else 
         Simpson := Simpson + 2*fvals[j] 
       fi 
     od; 
     RETURN(Simpson*h/3) 
     end:
 

Here are two examples:

> read `/home/fac/little/Classes/Numer-98/ListSimpson.map`; 
> fvals:=[1.4, 2.5, 2.3, 3.4, 5.6]; 
 
                  fvals := [1.4, 2.5, 2.3, 3.4, 5.6] 
With equally spaced x-values, h = .9: 
 
> ListSimpson(fvals,5,.9); 
 
                             10.56000000
 

Another example showing close agreement with the actual value of an integral – data points sampled from the sin function at \(x = .1\) to \(.9\) at intervals of \(.1\):

> fvals:=[seq(sin(.1+(i-1)*(.1)),i=1..9)]; 
 
  fvals := [.09983341665, .1986693308, .2955202067, .3894183423, 
 
        .4794255386, .5646424734, .6442176872, .7173560909, 
 
        .7833269096] 
 
> ListSimpson(fvals,9,.1); 
 
                             .3733944047 
 
> int(sin(x),x=0.1..0.9); 
 
                             .3733941970
 

For information on other similar methods, you can consult almost any textbook on numerical analysis. Hope this information is useful.

6.54.3 Michael Brennan (12.11.98)

All you need do is define a function, say \(f(x)\) using the experimental data. Then you can apply either the command simpson or trapezoid from the student package to find approximates to the area. Of course I am assuming the data you have is in the form of a table where the readings are taken repeatedly over a period of time.##e.g.

Suppose you are trying to find the area of a swimming pool where you are given the following measurements for the width of the pool, where the readings are taken every two metres.

Area of a Swimming Pool. First enter the data as a list.

> with(student); 
> L:=[0,6.2,7.2,6.8,5.6,5.0,4.8,4.8,0]: 
 
> for k from 0 to 8 do w(2*k):=L[k+1] od; 
 
> trapezoid(w(x),x=0..16,8); 
 
> evalf(%); 
 
> simpson(w(x),x=0..16,8); 
 
> evalf(%);
 

We would require about 85 m^2 of material to cover pool.

6.54.4 Helmut Kahovec (13.11.98)

If you have a mathematical model then fit it to the data. If the regression is linear then see the online help page for stats[fit,leastsquare]. If the regression is nonlinear then see the article

     Fit for Anything? 
     The Maple Reporter, Summer 1998 Edition 
     pp. 6-7
 

You may download The Maple Reporter from http://www.maplesoft.com. Afterwards transform the fit into a procedure using unapply() and (numerically) integrate this procedure using evalf/Int.

If you don’t have a mathematical model then I would suggest using splines to convert the data points into a (piecewise defined) procedure. See the online help page for spline(). If you have many points then let the degree of the segment polynomials be 1 (i.e., use lines), otherwise use cubic splines. However, if there are 30+ data points, the calculation of the splines may take some time, anyway. Afterwards (numerically) integrate the piecewise defined procedure within each interval separately using evalf/Int and add these areas using add(). Below I give you two examples of this.

> restart; 
First example: 
 
> f1:=x->exp(-x); 
 
                          f1 := x -> exp(-x) 
 
> evalf(Int( f1(x), x=0..2 )); 
 
                             .8646647168 
 
> l1:=[seq( [i/5,f1(i/5)], i=0..10 )]: 
> plot(l1,style=POINT); 
 
> readlib(spline): 
> f1s:=spline( 
>   map2(op,1,l1), 
>   map2(op,2,l1), 
>   x, 
>   3 
> ): 
> add( 
>   evalf(Int( f1s, x=l1[i][1]..l1[i+1][1] )), 
>   i=1..(nops(l1)-1) 
> ); 
 
                             .8648805898
 

Second example:

> f2:=x->ln(1+x); 
 
                         f2 := x -> ln(1 + x) 
 
> evalf(Int( f2(x), x=0..3 )); 
 
                             2.545177444 
 
> l2:=[seq([i/10,f2(i/10)],i=0..30)]: 
> plot(l2,style=POINT);
 

If you don’t need the splines explicitly you may define a procedure that calculates the area from this list of data points at once:

> AreaSpline:=proc(l::list(list),degr::integer) 
> local x,i; 
>   add( 
>     evalf( 
>       Int( 
>         readlib(spline)( 
>           map2(op,1,l), 
>           map2(op,2,l), 
>           x, 
>           degr 
>         ), 
>         x=l[i][1]..l[i+1][1] 
>       ) 
>     ), 
>     i=1..(nops(l)-1) 
>   ) 
> end: 
> AreaSpline(l2,1); 
 
                             2.544552717
 

6.54.5 Willard, Daniel (16.11.98)

Elementary calculus texts define integrals as the limit of a process of summing adjacent rectangles whose heights are the location of the left-hand ordinate (a data point), or the right-hand ordinate, or some average of the two. Simpson’s Rule. Trapezoidal interpolation.

Splines (cubic fits to successive sets of data points). See the math section of the CRC Handbook of Physics and Chemistry, or Abramowitz and Stegun Handbook of Mathematical Functions.