In
numerical analysis, numerical integration constitutes a broad family of
algorithms for calculating the numerical value of a definite integral, and by
extension, the term is also sometimes used to describe the numerical solution
of differential equations. The term numerical quadrature (often abbreviated to
quadrature) is more or less a synonym for numerical integration, especially as
applied to one-dimensional integrals.
The
basic problem in numerical integration is to compute an approximate solution to
a definite integral
to
a given degree of accuracy. If f(x) is a smooth function integrated over a
small number of dimensions, and the domain of integration is bounded, there are
many methods for approximating the integral to the desired precision.
An
important part of the analysis of any numerical integration method is to study
the behavior of the approximation error as a function of the number of
integrand evaluations. A method that yields a small error for a small number of
evaluations is usually considered superior. Reducing the number of evaluations
of the integrand reduces the number of arithmetic operations involved, and
therefore reduces the total round-off error. Also, each evaluation takes time,
and the integrand may be arbitrarily complicated.
A
'brute force' kind of numerical integration can be done, if the integrand is
reasonably well-behaved, by evaluating the integrand with very small
increments.
A
large class of quadrature rules can be derived by constructing interpolating
functions that are easy to integrate.
Typically
these interpolating functions are polynomials.
The
simplest method of this type is to let the interpolating function be a constant
function (a polynomial of degree zero) that passes through the point ((a+b)/2,
f((a+b)/2)). This is called the midpoint rule or rectangle rule.
The
interpolating function may be an affine function (a polynomial of degree 1)
that passes through the points (a, f(a)) and (b, f(b)). This is called the trapezoidal rule.
For
either one of these rules, we can make a more accurate approximation by
breaking up the interval [a, b] into some number n of subintervals, computing
an approximation for each subinterval, then adding up all the results. This is
called a composite rule, extended rule, or iterated rule. For
example, the composite trapezoidal rule can be stated as
where
the subintervals have the form [k h, (k+1) h], with h = (b−a)/n and k = 0, 1,
2, ..., n−1.
Interpolation
with polynomials evaluated at equally spaced points in [a, b] yields the
Newton–Cotes formulas, of which the rectangle rule and the trapezoidal rule are
examples. Simpson's rule is a Newton-Cotes formula for approximating the integral of a function 
using quadratic polynomials (i.e., parabolic arcs instead of the straight line segments used in the trapezoidal rule).
Quadrature
rules with equally spaced points have the very convenient property of
nesting. The corresponding rule with each interval subdivided includes all the
current points, so those integrand values can be re-used.
In sport and exercise biomechanics the data are most commonly available for wich the process of integration is appropiate are data from force platforms and accelerometers. The force data from force plataforms can be used to compute acceleration following Newton's second law (F = m * a). In this example, the process of integration enables the velocity to be obtained, and the process can be repeated in this velocity data in order to obtain displacement. This method can be applied in other type of biomechanics data, for example in EMG, like contractile impulse of activation. Integration of this type of data is best done using numerical integration.
The numerical integration method most commonly found in the literature in the area of biomechanics is the Trapezoidal. However, studies such as Won (2000) shows that we must have concern for the error produced by the choice of the integration method to be used and the method of Simpson is the one with the lowest error.
Simpson's rule can be derived by integrating a third-order Lagrange interpolating polynomial fit to the function at three equally spaced points. In particular, let the function 
be tabulated at points 
, 
, and 
equally spaced by distance 
, and denote 
. Then Simpson's rule states that
Since it uses quadratic polynomials to approximate functions, Simpson's rule actually gives exact results when approximating integrals of polynomials up to cubic degree. For example, consider 
(black curve) on the interval ![[0,pi/2]](https://lh3.googleusercontent.com/blogger_img_proxy/AEn0k_vnmiN-E1bdhnrhfsfZri01pp8fwKyq137l3fRT-SLhOBBD-xBrPnRNYrmFOfLWchc5CpPo9gJlFk7Y_me4suxn_CG7abUMYpUob1QfXDqjjDcuFYvtdO6PeQpY_yn9RoK7XR9fsssMoilh=s0-d)
, so that 
, 
, and 
.
Then Simpson's rule (which corresponds to the area under the blue curve obtained from the third-order interpolating polynomial) gives
whereas the trapezoidal rule (area under the red curve) gives 
and the actual answer is 1. In exact form,
where the remainder term can be written as
with 
being some value of 
in the interval ![[x_0,x_2]](https://lh3.googleusercontent.com/blogger_img_proxy/AEn0k_sUYGVbJmo0kSCQC5B3PSJ82DynOR_RxEZwn6KUD5QUCLJKZI0-IJ-R7GCPz4APZYRpFn2yIZ7crr7zVIG-WCy-RzCCaQToVOpDIPyP6EOKIQIj9D6EDLunIDRKQEpP5XlCWdMEqNjWjjgy=s0-d)
.
An extended version of the rule can be written for 
tabulated at 
, 
, ..., 
as
where the remainder term is
for some ![x^* in [x_0,x_(2n)]](https://lh3.googleusercontent.com/blogger_img_proxy/AEn0k_uGBs5lvNZVF7z3bXd6RuuqGLOG7FrzsWIj8ryspwBsVYGe6NvoSdrot2cAzGDDzUQRZkkJBFQoG9vXwvkTriH0xLKpFjnT3vFLmtoh8PIGeG-H4Rj69yM6v7GIrP052NL7pOUyQLmGiEMt=s0-d)
.
Here's a link to download Matlab Code with a sub-program implemented in environment Matlab© to get the integration value of your signal with Simpson's Rule.
To
use the function provided just save the file "simpson_gbiomech.m" in the
same folder as the data to be analyzed. Then just call in your Matlab
routine sub-program as follows in the example below:
Z=simpson_gbiomech(Y); - computes an approximation of the integral of Y via the Simpson's method (with unit spacing). To compute the integral for spacing different from one, multiply Z by the spacing increment.
Z=simpson_gbiomech(X,Y); - computes the integral of Y with respect to X using the Simpson's rule. X and Y must be vectors of the same length, or X must be a column vector and Y an array whose first non-singleton dimension is length(X). SIMPS operates along this dimension.
Note: you can specify the window which will be modified by inserting an index on the input variables.
References:
- Abramowitz, M. and Stegun, I. A. (Eds.). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. New York: Dover, p. 886, 1972.
- Davis, P. J.; Rabinowitz, P.. Methods of Numerical Integration.
- Horwitz, A. "A Version of Simpson's Rule for Multiple Integrals." J. Comput. Appl. Math. 134, 1-11, 2001.
- Jeffreys, H. and Jeffreys, B. S. Methods of Mathematical Physics, 3rd ed. Cambridge, England: Cambridge University Press, p. 286, 1988.
- Stoer J.; Bulirsch R.. Introduction to Numerical Analysis. New York: Springer-Verlag, 1980. (See Chapter 3.)
- Trott, M. The Mathematica GuideBook for Programming. New York: Springer-Verlag, p. 105, 2004. http://www.mathematicaguidebooks.org/.
- Whittaker, E. T. and Robinson, G. "The Trapezoidal and Parabolic Rules." The Calculus of Observations: A Treatise on Numerical Mathematics, 4th ed. New York: Dover, pp. 156-158, 1967.
- Won, T. The comparison of among different numerical integration methods on the biomechanics of sport countermovement jump. 18º International Symposium on Biomechanics in Sports. Hong Kong, China: 2000. Available: <https://ojs.ub.uni-konstanz.de/cpa/article/view/2270/2123>
Big Hug!
