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: michaelvdk@hotmail.nospam.com
(omit the nospam for the true address).
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.
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.
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
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.