|
|
| (620 intermediate revisions by more than 100 users not shown) |
| Line 1: |
Line 1: |
| [[Image:Integral as region under curve.svg|thumb|right|Numerical integration consists of finding numerical approximations for the value <math>S</math>]]
| | This is a preview for the new '''MathML rendering mode''' (with SVG fallback), which is availble in production for registered users. |
| 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 ordinary differential equations|numerical solution of differential equations]]. This article focuses on calculation of definite integrals. The term '''numerical quadrature''' (often abbreviated to [[Quadrature (mathematics)|''quadrature'']]) is more or less a synonym for ''numerical integration'', especially as applied to one-dimensional integrals. Numerical integration over more than one dimension is sometimes described as '''cubature''',<ref>{{MathWorld | urlname=Cubature | title=Cubature }}</ref> although the meaning of ''quadrature'' is understood for higher dimensional integration as well.
| |
|
| |
|
| The basic problem considered by numerical integration is to compute an approximate solution to a definite integral:
| | If you would like use the '''MathML''' rendering mode, you need a wikipedia user account that can be registered here [[https://en.wikipedia.org/wiki/Special:UserLogin/signup]] |
| | * Only registered users will be able to execute this rendering mode. |
| | * Note: you need not enter a email address (nor any other private information). Please do not use a password that you use elsewhere. |
|
| |
|
| :<math>\int_a^b\! f(x)\, dx.</math> | | Registered users will be able to choose between the following three rendering modes: |
|
| |
|
| If {{math|''f(x)''}} is a smooth well-behaved function, integrated over a small number of dimensions and the limits of integration are bounded, there are many methods of approximating the integral with arbitrary precision.
| | '''MathML''' |
| | :<math forcemathmode="mathml">E=mc^2</math> |
|
| |
|
| ==Reasons for numerical integration== | | <!--'''PNG''' (currently default in production) |
| | :<math forcemathmode="png">E=mc^2</math> |
|
| |
|
| There are several reasons for carrying out numerical integration.
| | '''source''' |
| The integrand ''f(x)'' may be known only at certain points,
| | :<math forcemathmode="source">E=mc^2</math> --> |
| such as obtained by [[sampling (statistics)|sampling]].
| |
| Some [[embedded systems]] and other computer applications may need numerical integration for this reason.
| |
|
| |
|
| A formula for the integrand may be known, but it may be difficult or impossible to find an [[antiderivative]] which is an [[elementary function]]. An example of such an integrand is ''f(x)'' = exp(−''x''<sup>''2''</sup>), the antiderivative of which (the [[error function]], times a constant) cannot be written in [[elementary form]].
| | <span style="color: red">Follow this [https://en.wikipedia.org/wiki/Special:Preferences#mw-prefsection-rendering link] to change your Math rendering settings.</span> You can also add a [https://en.wikipedia.org/wiki/Special:Preferences#mw-prefsection-rendering-skin Custom CSS] to force the MathML/SVG rendering or select different font families. See [https://www.mediawiki.org/wiki/Extension:Math#CSS_for_the_MathML_with_SVG_fallback_mode these examples]. |
|
| |
|
| It may be possible to find an antiderivative symbolically, but it may be easier to compute a numerical approximation than to compute the antiderivative. That may be the case if the antiderivative is given as an infinite series or product, or if its evaluation requires a [[special function]] which is not available.
| | ==Demos== |
|
| |
|
| ==Methods for one-dimensional integrals== | | Here are some [https://commons.wikimedia.org/w/index.php?title=Special:ListFiles/Frederic.wang demos]: |
|
| |
|
| Numerical integration methods can generally be described as combining evaluations of the integrand to get an approximation to the integral. The integrand is evaluated at a finite set of points called '''integration points''' and a weighted sum of these values is used to approximate the integral. The integration points and weights depend on the specific method used and the accuracy required from the approximation.
| |
|
| |
|
| 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.
| | * accessibility: |
| A method which yields a small error for a small number of evaluations is usually considered superior.
| | ** Safari + VoiceOver: [https://commons.wikimedia.org/wiki/File:VoiceOver-Mac-Safari.ogv video only], [[File:Voiceover-mathml-example-1.wav|thumb|Voiceover-mathml-example-1]], [[File:Voiceover-mathml-example-2.wav|thumb|Voiceover-mathml-example-2]], [[File:Voiceover-mathml-example-3.wav|thumb|Voiceover-mathml-example-3]], [[File:Voiceover-mathml-example-4.wav|thumb|Voiceover-mathml-example-4]], [[File:Voiceover-mathml-example-5.wav|thumb|Voiceover-mathml-example-5]], [[File:Voiceover-mathml-example-6.wav|thumb|Voiceover-mathml-example-6]], [[File:Voiceover-mathml-example-7.wav|thumb|Voiceover-mathml-example-7]] |
| Reducing the number of evaluations of the integrand reduces the number of arithmetic operations involved,
| | ** [https://commons.wikimedia.org/wiki/File:MathPlayer-Audio-Windows7-InternetExplorer.ogg Internet Explorer + MathPlayer (audio)] |
| and therefore reduces the total [[round-off error]].
| | ** [https://commons.wikimedia.org/wiki/File:MathPlayer-SynchronizedHighlighting-WIndows7-InternetExplorer.png Internet Explorer + MathPlayer (synchronized highlighting)] |
| Also,
| | ** [https://commons.wikimedia.org/wiki/File:MathPlayer-Braille-Windows7-InternetExplorer.png Internet Explorer + MathPlayer (braille)] |
| each evaluation takes time, and the integrand may be arbitrarily complicated.
| | ** NVDA+MathPlayer: [[File:Nvda-mathml-example-1.wav|thumb|Nvda-mathml-example-1]], [[File:Nvda-mathml-example-2.wav|thumb|Nvda-mathml-example-2]], [[File:Nvda-mathml-example-3.wav|thumb|Nvda-mathml-example-3]], [[File:Nvda-mathml-example-4.wav|thumb|Nvda-mathml-example-4]], [[File:Nvda-mathml-example-5.wav|thumb|Nvda-mathml-example-5]], [[File:Nvda-mathml-example-6.wav|thumb|Nvda-mathml-example-6]], [[File:Nvda-mathml-example-7.wav|thumb|Nvda-mathml-example-7]]. |
| | ** Orca: There is ongoing work, but no support at all at the moment [[File:Orca-mathml-example-1.wav|thumb|Orca-mathml-example-1]], [[File:Orca-mathml-example-2.wav|thumb|Orca-mathml-example-2]], [[File:Orca-mathml-example-3.wav|thumb|Orca-mathml-example-3]], [[File:Orca-mathml-example-4.wav|thumb|Orca-mathml-example-4]], [[File:Orca-mathml-example-5.wav|thumb|Orca-mathml-example-5]], [[File:Orca-mathml-example-6.wav|thumb|Orca-mathml-example-6]], [[File:Orca-mathml-example-7.wav|thumb|Orca-mathml-example-7]]. |
| | ** From our testing, ChromeVox and JAWS are not able to read the formulas generated by the MathML mode. |
|
| |
|
| A 'brute force' kind of numerical integration can be done, if the integrand is reasonably well-behaved (i.e. [[piecewise]] [[continuous function|continuous]] and of [[bounded variation]]), by evaluating the integrand with very small increments.
| | ==Test pages == |
|
| |
|
| === Quadrature rules based on interpolating functions ===
| | To test the '''MathML''', '''PNG''', and '''source''' rendering modes, please go to one of the following test pages: |
| | *[[Displaystyle]] |
| | *[[MathAxisAlignment]] |
| | *[[Styling]] |
| | *[[Linebreaking]] |
| | *[[Unique Ids]] |
| | *[[Help:Formula]] |
|
| |
|
| A large class of quadrature rules can be derived by constructing [[interpolation|interpolating]] functions which are easy to integrate. Typically these interpolating functions are [[polynomial]]s.
| | *[[Inputtypes|Inputtypes (private Wikis only)]] |
| | | *[[Url2Image|Url2Image (private Wikis only)]] |
| [[Image:Integration rectangle.svg|right|frame|Illustration of the rectangle rule.]] | | ==Bug reporting== |
| The simplest method of this type is to let the interpolating function be a constant function (a polynomial of degree zero) which passes through the point ((''a''+''b'')/2, ''f''((''a''+''b'')/2)). This is called the ''midpoint rule'' or ''[[rectangle method|rectangle rule]]''.
| | If you find any bugs, please report them at [https://bugzilla.wikimedia.org/enter_bug.cgi?product=MediaWiki%20extensions&component=Math&version=master&short_desc=Math-preview%20rendering%20problem Bugzilla], or write an email to math_bugs (at) ckurs (dot) de . |
| | |
| :<math>\int_a^b f(x)\,dx \approx (b-a) \, f\left(\frac{a+b}{2}\right).</math>
| |
| | |
| [[Image:Integration trapezoid.svg|right|frame|Illustration of the trapezoidal rule. ]] | |
| The interpolating function may be an [[affine function]] (a polynomial of degree 1)
| |
| which passes through the points (''a'', ''f''(''a'')) and (''b'', ''f''(''b'')).
| |
| This is called the ''[[trapezoidal rule]]''.
| |
| | |
| :<math>\int_a^b f(x)\,dx \approx (b-a) \, \frac{f(a) + f(b)}{2}.</math>
| |
| | |
| [[Image:Integration simpson.png|right|frame|Illustration of Simpson's 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
| |
| | |
| :<math>\int_a^b f(x)\,dx \approx \frac{b-a}{n} \left( {f(a) \over 2} + \sum_{k=1}^{n-1} f \left( a+k \frac{b-a}{n} \right) + {f(b) \over 2} \right)</math>
| |
| | |
| 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]], which is based on a polynomial of order 2, is also a Newton–Cotes formula.
| |
| | |
| Quadrature rules with equally spaced points have the very convenient property of <em>nesting</em>. The corresponding rule with each interval subdivided includes all the current points, so those integrand values can be re-used.
| |
| | |
| If we allow the intervals between interpolation points to vary, we find another group of quadrature formulas, such as the [[Gaussian quadrature]] formulas. A Gaussian quadrature rule is typically more accurate than a Newton–Cotes rule which requires the same number of function evaluations, if the integrand is [[Smooth function|smooth]] (i.e., if it is sufficiently differentiable). Other quadrature methods with varying intervals include [[Clenshaw–Curtis quadrature]] (also called Fejér quadrature) methods, which do nest. | |
| | |
| Gaussian quadrature rules do not nest, but the related [[Gauss–Kronrod quadrature formula]]s do.
| |
| | |
| === Adaptive algorithms ===
| |
| {{details|Adaptive quadrature}}
| |
| | |
| If ''f(x)'' does not have many derivatives at all points, or if the derivatives become large, then Gaussian quadrature is often insufficient. In this case, an algorithm similar to the following will perform better:
| |
| | |
| | |
| <source lang=python>
| |
| def calculate_definite_integral_of_f(f, initial_step_size):
| |
| '''
| |
| This algorithm calculates the definite integral of a function
| |
| from 0 to 1, adaptively, by choosing smaller steps near
| |
| problematic points.
| |
| '''
| |
| x = 0.0
| |
| h = initial_step_size
| |
| accumulator = 0.0
| |
| while x < 1.0:
| |
| if x + h > 1.0:
| |
| h = 1.0 - x
| |
| quad_this_step =
| |
| if error_too_big_in_quadrature_of_over_range(f, [x,x+h]):
| |
| h = make_h_smaller(h)
| |
| else:
| |
| accumulator += quadrature_of_f_over_range(f, [x,x+h])
| |
| x += h
| |
| if error_too_small_in_quadrature_of_over_range(f, [x,x+h]):
| |
| h = make_h_larger(h) # Avoid wasting time on tiny steps.
| |
| return accumulator
| |
| </source>
| |
| | |
| Some details of the algorithm require careful thought. For many cases, estimating the error from quadrature over an interval for a function ''f''(''x'') isn't obvious. One popular solution is to use two different rules of quadrature, and use their difference as an estimate of the error from quadrature. The other problem is deciding what "too large" or "very small" signify. A <em>local</em> criterion for "too large" is that the quadrature error should not be larger than ''t'' · ''h'' where ''t'', a real number, is the tolerance we wish to set for global error. Then again, if ''h'' is already tiny, it may not be worthwhile to make it even smaller even if the quadrature error is apparently large. A <em>global</em> criterion is that the sum of errors on all the intervals should be less than ''t''. This type of error analysis is usually called "a posteriori" since we compute the error after having computed the approximation.
| |
| | |
| Heuristics for adaptive quadrature are discussed by Forsythe et al. (Section 5.4).
| |
| | |
| === Extrapolation methods ===
| |
| | |
| The accuracy of a quadrature rule of the Newton-Cotes type is generally a function of the number of evaluation points.
| |
| The result is usually more accurate as number of evaluation points increases,
| |
| or, equivalently, as the width of the step size between the points decreases.
| |
| It is natural to ask what the result would be if the step size were allowed to approach zero.
| |
| This can be answered by extrapolating the result from two or more nonzero step sizes, using [[series acceleration]] methods such as [[Richardson extrapolation]].
| |
| The extrapolation function may be a [[polynomial]] or [[rational function]].
| |
| Extrapolation methods are described in more detail by Stoer and Bulirsch (Section 3.4) and are implemented in many of the routines in the [[QUADPACK]] library.
| |
| | |
| === Conservative (a priori) error estimation ===
| |
| | |
| Let ''f'' have a bounded first derivative over [''a'',''b'']. The [[mean value theorem]] for ''f'', where ''x'' < ''b'', gives
| |
| | |
| : <math>(x - a) f'(y_x) = f(x) - f(a)\,</math>
| |
| | |
| for some ''y<sub>x</sub>'' in [''a'',''x''] depending on ''x''. If we integrate in ''x'' from ''a'' to ''b'' on both sides and take the absolute values, we obtain
| |
| | |
| : <math>\left| \int_a^b f(x)\,dx - (b - a) f(a) \right|
| |
| = \left| \int_a^b (x - a) f'(y_x)\, dx \right|</math>
| |
| | |
| We can further approximate the integral on the right-hand side by bringing the absolute value into the integrand, and replacing the term in ''f' '' by an upper bound:
| |
| | |
| : <math>\left| \int_a^b f(x)\,dx - (b - a) f(a) \right| \leq {(b - a)^2 \over 2} \sup_{a \leq x \leq b} \left| f'(x) \right|</math> (**)
| |
| | |
| (See [[supremum]].) Hence, if we approximate the integral ∫<sub>''a''</sub><sup>''b''</sup> ''f''(''x'') d''x'' by the quadrature rule (''b'' − ''a'')''f''(''a'') our error is no greater than the right hand side of (**). We can convert this into an error analysis for the Riemann sum (*), giving an upper bound of
| |
| | |
| : <math>{n^{-1} \over 2} \sup_{0 \leq x \leq 1} \left| f'(x) \right|</math>
| |
| | |
| for the [[error term]] of that particular approximation. (Note that this is precisely the error we calculated for the example <math>f(x) = x</math>.) Using more derivatives, and by tweaking the quadrature, we can do a similar error analysis using a [[Taylor series]] (using a partial sum with remainder term) for ''f''. This error analysis gives a strict upper bound on the error, if the derivatives of ''f'' are available.
| |
| | |
| This integration method can be combined with [[interval arithmetic]] to produce [[computer proof]]s and ''verified'' calculations.
| |
| | |
| === Integrals over infinite intervals ===
| |
| | |
| ==== Infinite intervals ====
| |
| One way to calculate an integral over infinite interval,
| |
| | |
| :<math>
| |
| \int_{-\infty}^{+\infty}f(x) \, dx,
| |
| </math>
| |
| | |
| is to transform it into an integral over a finite interval by any one of several possible changes of variables, for example:
| |
| | |
| :<math>
| |
| \int_{-\infty}^{+\infty} f(x) \, dx = \int_{-1}^{+1} f\left( \frac{t}{1-t^2} \right) \frac{1+t^2}{(1-t^2)^2} \, dt,
| |
| </math>
| |
| | |
| The integral over finite interval can then be evaluated by ordinary integration methods.
| |
| | |
| ==== Half-infinite intervals ====
| |
| An integral over a half-infinite interval can likewise be transformed into an integral over a finite interval by any one of several possible changes of variables, for example:
| |
| | |
| :<math>
| |
| \int_a^{+\infty}f(x) \, dx =\int_0^1 f\left(a + \frac{1-t}{t}\right) \frac{dt}{t^2} .</math>
| |
| | |
| Similarly,
| |
| | |
| :<math>
| |
| \int_{-\infty}^a f(x) \, dx = \int_0^1 f\left(a - \frac{1-t}{t}\right) \frac{dt}{t^2}</math>
| |
| | |
| == Multidimensional integrals ==
| |
| | |
| The quadrature rules discussed so far are all designed to compute one-dimensional integrals.
| |
| To compute integrals in multiple dimensions,
| |
| one approach is to phrase the multiple integral as repeated one-dimensional integrals by appealing to [[Fubini's theorem]].
| |
| This approach requires the function evaluations to [[exponential growth|grow exponentially]] as the number of dimensions increases. Two methods are known to overcome this so-called ''[[curse of dimensionality]]''.
| |
| | |
| === Monte Carlo ===
| |
| | |
| {{main|Monte Carlo integration}}
| |
| | |
| [[Monte Carlo method]]s and [[quasi-Monte Carlo method]]s are easy to apply to multi-dimensional integrals,
| |
| and may yield greater accuracy for the same number of function evaluations than repeated integrations using one-dimensional methods.
| |
| | |
| A large class of useful Monte Carlo methods are the so-called [[Markov chain Monte Carlo]] algorithms,
| |
| which include the [[Metropolis-Hastings algorithm]] and [[Gibbs sampling]].
| |
| | |
| === Sparse grids ===
| |
| [[Sparse grid]]s were originally developed by Smolyak for the quadrature of high dimensional functions. The method is always based on a one dimensional quadrature rule, but performs a more sophisticated combination of univariate results.
| |
| | |
| == Connection with differential equations ==
| |
| | |
| The problem of evaluating the integral
| |
| :<math>\int_a^b f(x)\, dx</math>
| |
| can be reduced to an [[initial value problem]] for an [[ordinary differential equation]]. If the above integral is denoted by ''I''(''b''), then the function ''I'' satisfies
| |
| :<math> I'(x) = f(x), \quad I(a) = 0. </math>
| |
| Methods developed for ordinary differential equations, such as [[Runge–Kutta methods]], can be applied to the restated problem and thus be used to evaluate the integral. For instance, the standard fourth-order Runge–Kutta method applied to the differential equation yields Simpson's rule from above.
| |
| | |
| The differential equation ''I'' ' (''x'') = ƒ(''x'') has a special form: the right-hand side contains only the dependent variable (here ''x'') and not the independent variable (here ''I''). This simplifies the theory and algorithms considerably. The problem of evaluating integrals is thus best studied in its own right.
| |
| | |
| ==See also==
| |
| * [[Numerical ordinary differential equations]]
| |
| * [[Truncation error (numerical integration)]]
| |
| * [[Clenshaw–Curtis quadrature]]
| |
| * [[Gauss-Kronrod quadrature]]
| |
| * [[Riemann Sum]] or [[Riemann Integral]]
| |
| * [[Trapezoidal Rule]]
| |
| | |
| == References ==
| |
| {{reflist}}
| |
| * [[Philip J. Davis]] and [[Philip Rabinowitz (mathematician)|Philip Rabinowitz]], ''Methods of Numerical Integration''.
| |
| * George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler. ''Computer Methods for Mathematical Computations''. Englewood Cliffs, NJ: Prentice-Hall, 1977. ''(See Chapter 5.)''
| |
| * {{Citation |last1=Press|first1=WH|last2=Teukolsky|first2=SA|last3=Vetterling|first3=WT|last4=Flannery|first4=BP|year=2007|title=Numerical Recipes: The Art of Scientific Computing|edition=3rd|publisher=Cambridge University Press| publication-place=New York|isbn=978-0-521-88068-8|chapter=Chapter 4. Integration of Functions|chapter-url=http://apps.nrbook.com/empanel/index.html?pg=155}}
| |
| * Josef Stoer and Roland Bulirsch. ''Introduction to Numerical Analysis''. New York: Springer-Verlag, 1980. ''(See Chapter 3.)''
| |
| | |
| ==External links==
| |
| * [http://numericalmethods.eng.usf.edu/mws/gen/07int/index.html Integration: Background, Simulations, etc.] at Holistic Numerical Methods Institute
| |
| | |
| === Free software for numerical integration ===
| |
| | |
| Numerical integration is one of the most intensively studied problems in numerical analysis.
| |
| Of the many software implementations, we list a few [[free and open source software]] packages here:
| |
| | |
| * [[QUADPACK]] (part of SLATEC): description [http://www.netlib.org/slatec/src/qpdoc.f], source code [http://www.netlib.org/slatec/src]. QUADPACK is a collection of algorithms, in Fortran, for numerical integration based on Gaussian quadrature.
| |
| * [http://openopt.org/interalg interalg]: a solver from [[OpenOpt]]/[[FuncDesigner]] frameworks, based on interval analysis, '''guaranteed precision''', license: BSD (free for any purposes)
| |
| * [http://www.gnu.org/software/gsl/ GSL]: The GNU Scientific Library (GSL) is a numerical library written in C which provides a wide range of mathematical routines, like Monte Carlo integration.
| |
| * Numerical integration algorithms are found in [[Guide to Available Mathematical Software|GAMS]] class [http://gams.nist.gov/serve.cgi/Class/H2 H2].
| |
| * [http://www.alglib.net/integral/ ALGLIB] is a collection of algorithms, in C# / C++ / Delphi / Visual Basic / etc., for numerical integration (includes Bulirsch-Stoer and Runge-Kutta integrators).
| |
| * [http://www.feynarts.de/cuba/ Cuba] is a free-software library of several multi-dimensional integration algorithms.
| |
| * [http://ab-initio.mit.edu/wiki/index.php/Cubature Cubature] code for adaptive multi-dimensional integration.
| |
| * [http://www.scilab.org/ Scilab] is an open source software under GPL license, providing powerful features including numerical integration.
| |
| | |
| [[Category:Numerical analysis]]
| |
| [[Category:Numerical integration (quadrature)|*]]
| |
| [[Category:Articles with example Python code]]
| |
| | |
| [[ar:تكامل عددي]]
| |
| [[ca:Integració numèrica]]
| |
| [[de:Numerische Integration]]
| |
| [[es:Integración numérica]]
| |
| [[fr:Calcul numérique d'une intégrale]]
| |
| [[ko:수치적분]]
| |
| [[it:Integrazione numerica]]
| |
| [[kk:Сандық интегралдау]]
| |
| [[he:שיטות נומריות לחישוב אינטגרלים מסוימים]]
| |
| [[hu:Numerikus integrálás]]
| |
| [[nl:Numerieke integratie]]
| |
| [[ja:数値積分]]
| |
| [[pl:Całkowanie numeryczne]]
| |
| [[pt:Integração numérica]]
| |
| [[ru:Численное интегрирование]]
| |
| [[sr:Нумеричка интеграција]]
| |
| [[sv:Numerisk integrering]]
| |
| [[uk:Чисельне інтегрування]]
| |
| [[zh:數值積分]]
| |
