--- title: "NIntegrate" language: "en" type: "Symbol" summary: "NIntegrate[f, {x, xmin, xmax}] gives a numerical approximation to the integral \\[Integral]_xmin^xmax\\ f\\ d x. NIntegrate[f, {x, xmin, xmax}, {y, ymin, ymax}, ...] gives a numerical approximation to the multiple integral \\[Integral]_xmin^xmaxd x \\[Integral]_ymin^ymaxd y\\ \\ ...\\ f. NIntegrate[f, {x, y, ...} \\[Element] reg] integrates over the geometric region reg." keywords: - adaptive Monte Carlo - adaptive Monte Carlo integration - adaptive quasi Monte Carlo - adaptive quasi Monte Carlo integration - product cubature - rule - Clenshaw–Curtis integration - Clenshaw–Curtis rule - composite quadrature - cubature - double-exponential - double exponential integration - Gaussian quadrature - Gauss integration - Gauss–Kronrod integration - Gauss–Kronrod rule - Genz–Malik algorithm - Genz–Malik cubature - global adaptive - global adaptive integration - integration - Kronrod points - Las Vegas integration - Lobatto–Kronrod integration - Lobatto–Kronrod rule - local adaptive - local adaptive integration - Monte Carlo - Monte Carlo integration - multidimensional rule - multi-panel quadrature - multipanel rule - Newton–Cotes integration - Newton–Cotes rule - nint - numerical integration - quadrature - quasi Monte Carlo - quasi Monte Carlo integration - symbolic processing - trapezoidal integration - trapezoidal rule - INT_2D - INT_3D - QROMB - QROMO - QSIMP - dblquad - quad - quadl - quadv - trapz - triplequad canonical_url: "https://reference.wolfram.com/language/ref/NIntegrate.html" source: "Wolfram Language Documentation" related_guides: - title: "Calculus" link: "https://reference.wolfram.com/language/guide/Calculus.en.md" - title: "Geometric Computation" link: "https://reference.wolfram.com/language/guide/GeometricComputation.en.md" - title: "Functions of Complex Variables" link: "https://reference.wolfram.com/language/guide/FunctionsOfComplexVariables.en.md" - title: "Solvers over Regions" link: "https://reference.wolfram.com/language/guide/GeometricSolvers.en.md" - title: "Symbolic Vectors, Matrices and Arrays" link: "https://reference.wolfram.com/language/guide/SymbolicArrays.en.md" - title: "Numerical Evaluation & Precision" link: "https://reference.wolfram.com/language/guide/NumericalEvaluationAndPrecision.en.md" - title: "Mesh-Based Geometric Regions" link: "https://reference.wolfram.com/language/guide/MeshRegions.en.md" - title: "Polygons" link: "https://reference.wolfram.com/language/guide/Polygons.en.md" - title: "Polyhedra" link: "https://reference.wolfram.com/language/guide/Polyhedra.en.md" related_functions: - title: "NDSolve" link: "https://reference.wolfram.com/language/ref/NDSolve.en.md" - title: "NSum" link: "https://reference.wolfram.com/language/ref/NSum.en.md" - title: "Integrate" link: "https://reference.wolfram.com/language/ref/Integrate.en.md" - title: "AsymptoticIntegrate" link: "https://reference.wolfram.com/language/ref/AsymptoticIntegrate.en.md" - title: "NLineIntegrate" link: "https://reference.wolfram.com/language/ref/NLineIntegrate.en.md" - title: "NSurfaceIntegrate" link: "https://reference.wolfram.com/language/ref/NSurfaceIntegrate.en.md" - title: "NContourIntegrate" link: "https://reference.wolfram.com/language/ref/NContourIntegrate.en.md" - title: "NExpectation" link: "https://reference.wolfram.com/language/ref/NExpectation.en.md" - title: "NProbability" link: "https://reference.wolfram.com/language/ref/NProbability.en.md" - title: "ArcLength" link: "https://reference.wolfram.com/language/ref/ArcLength.en.md" - title: "Area" link: "https://reference.wolfram.com/language/ref/Area.en.md" - title: "Volume" link: "https://reference.wolfram.com/language/ref/Volume.en.md" - title: "MomentOfInertia" link: "https://reference.wolfram.com/language/ref/MomentOfInertia.en.md" related_tutorials: - title: "Advanced Numerical Integration in the Wolfram Language" link: "https://reference.wolfram.com/language/tutorial/NIntegrateOverview.en.md" - title: "Numerical Mathematics: Basic Operations" link: "https://reference.wolfram.com/language/tutorial/NumericalCalculations.en.md" - title: "Numerical Integration" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#12104" - title: "Numerical Sums, Products, and Integrals" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#3848" - title: "Numerical Mathematics in the Wolfram Language" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#20411" - title: "The Uncertainties of Numerical Mathematics" link: "https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#14523" - title: "Implementation notes: Numerical and Related Functions" link: "https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#23656" --- # NIntegrate NIntegrate[f, {x, xmin, xmax}] gives a numerical approximation to the integral $\int _{x_{\text{\textit{min}}}}^{x_{\text{\textit{max}}}} f dx$. NIntegrate[f, {x, xmin, xmax}, {y, ymin, ymax}, …] gives a numerical approximation to the multiple integral $\int _{x_{\text{\textit{min}}}}^{x_{\text{\textit{max}}}}dx\int _{y_{\text{\textit{min}}}}^{y_{\text{\textit{max}}}}dy \ldots f$. NIntegrate[f, {x, y, …}∈reg] integrates over the geometric region reg. ## Details and Options * Multiple integrals use a variant of the standard iterator notation. The first variable given corresponds to the outermost integral and is done last. * ``NIntegrate`` by default tests for singularities at the boundaries of the integration region and at the boundaries of regions specified by settings for the ``Exclusions`` option. * ``NIntegrate[f, {x, x0, x1, …, xk}]`` tests for singularities in a one-dimensional integral at each of the intermediate points ``xi``. If there are no singularities, the result is equivalent to an integral from ``x0`` to ``xk``. You can use complex numbers ``xi`` to specify an integration contour in the complex plane. * The following options can be given: | | | | | :----------------- | :--------------- | :------------------------------------------------- | | AccuracyGoal | Infinity | digits of absolute accuracy sought | | EvaluationMonitor | None | expression to evaluate whenever expr is evaluated | | Exclusions | None | parts of the integration region to exclude | | MaxPoints | Automatic | maximum total number of sample points | | MaxRecursion | Automatic | maximum number of recursive subdivisions | | Method | Automatic | method to use | | MinRecursion | 0 | minimum number of recursive subdivisions | | PrecisionGoal | Automatic | digits of precision sought | | WorkingPrecision | MachinePrecision | the precision used in internal computations | * ``NIntegrate`` usually uses adaptive algorithms, which recursively subdivide the integration region as needed. ``MinRecursion`` specifies the minimum number of recursive subdivisions to try. ``MaxRecursion`` gives the maximum number. * ``NIntegrate`` usually continues doing subdivisions until the error estimate it gets implies that the final result achieves either the ``AccuracyGoal`` or the ``PrecisionGoal`` specified. * For low-dimensional integrals, the default setting for ``PrecisionGoal`` is related to ``WorkingPrecision``. For high-dimensional integrals, it is typically taken to be a fixed value, usually ``2``. * You should realize that with sufficiently pathological functions, the algorithms used by ``NIntegrate`` can give wrong answers. In most cases, you can test the answer by looking at its sensitivity to changes in the setting of options for ``NIntegrate``. * Possible explicit settings for the ``Method`` option include: | | | | :------------------------ | :------------------------------------- | | "GlobalAdaptive" | global adaptive integration strategy | | "LocalAdaptive" | local adaptive integration strategy | | "DoubleExponential" | double exponential quadrature | | "MonteCarlo" | Monte Carlo integration | | "AdaptiveMonteCarlo" | adaptive Monte Carlo integration | | "QuasiMonteCarlo" | quasi Monte Carlo integration | | "AdaptiveQuasiMonteCarlo" | adaptive quasi Monte Carlo integration | * With ``Method -> {"strategy", Method -> "rule"}`` or ``Method -> {"strategy", Method -> {rule1, rule2, …}}`` possible strategy methods include: | | | | :--------------- | :-------------------------------------------- | | "GlobalAdaptive" | subdivide based on global error estimates | | "LocalAdaptive" | subdivide based only on local error estimates | * Methods used as rules include: | | | | :--------------------- | :------------------------------------------ | | "CartesianRule" | multidimensional product of rules | | "ClenshawCurtisRule" | Clenshaw–Curtis rule | | "GaussKronrodRule" | Gauss points with Kronrod extension | | "LevinRule" | Levin-type oscillatory rules | | "LobattoKronrodRule" | Gauss–Lobatto points with Kronrod extension | | "MultidimensionalRule" | multidimensional symmetric rule | | "MultipanelRule" | combination of 1D rules | | "NewtonCotesRule" | Newton–Cotes rule | | "RiemannRule" | Riemann sum rule | | "TrapezoidalRule" | uniform points in one dimension | * With the setting ``Method -> "rule"``, the strategy method will be selected automatically. * Additional method suboptions can be given in the form ``Method -> {…, opts}``. * ``NIntegrate`` symbolically analyzes its input to transform oscillatory and other integrands, subdivide piecewise functions, and select optimal algorithms. * The method suboption ``"SymbolicProcessing"`` specifies the maximum number of seconds for which to attempt performing symbolic analysis of the integrand. * ``N[Integrate[…]]`` calls ``NIntegrate`` for integrals that cannot be done symbolically. * ``NIntegrate`` first localizes the values of all variables, then evaluates ``f`` with the variables being symbolic, and then repeatedly evaluates the result numerically. * ``NIntegrate`` has attribute ``HoldAll`` and effectively uses ``Block`` to localize variables. --- ## Examples (134) ### Basic Examples (5) Compute a numerical integral: ```wl In[1]:= NIntegrate[Sin[Sin[x]], {x, 0, 2}] Out[1]= 1.24706 ``` --- Compute a multidimensional integral (with singularity at the origin): ```wl In[1]:= NIntegrate[(1/Sqrt[x^1 + y^2 + z^3]), {x, 0, 1}, {y, 0, 1}, {z, 0, 1}] Out[1]= 1.0885 ``` --- Compute areas and volumes for implicitly defined regions: ```wl In[1]:= NIntegrate[Boole[1 / 4 ≤ x ^ 2 + (2y) ^ 2 ≤ 1∧y > x], {x, -1, 1}, {y, -1, 1}] Out[1]= 0.589049 In[2]:= RegionPlot[1 / 4 ≤ x ^ 2 + (2y) ^ 2 ≤ 1∧y > x, {x, -1, 1}, {y, -1, 1}] Out[2]= [image] ``` --- Integrate over any region: ```wl In[1]:= ℛ = ConvexHullMesh[RandomReal[1, {50, 3}]] Out[1]= [image] In[2]:= NIntegrate[x ^ 2y ^ 2z ^ 2, {x, y, z}∈ℛ] Out[2]= 0.0124094 ``` --- Integrate oscillatory and other complicated functions: ```wl In[1]:= NIntegrate[(1/x)Cos[Log[x] / x], {x, 0, 1}] Out[1]= 0.323367 In[2]:= Plot[(1/x)Cos[Log[x] / x], {x, 0, 1}, PlotRange -> 15] Out[2]= [image] ``` ### Scope (56) #### Basic Uses (12) Integrate functions over a finite real range: ```wl In[1]:= NIntegrate[Exp[Exp[-x]], {x, 0, 5}] Out[1]= 6.31115 In[2]:= NIntegrate[x ^ x, {x, 1, 5}] Out[2]= 1241.03 ``` --- An infinite real range: ```wl In[1]:= NIntegrate[(1/1 + x ^ 2), {x, 0, ∞}] Out[1]= 1.5708 In[2]:= NIntegrate[Exp[I x ^ 2], {x, -∞, ∞}] Out[2]= 1.25331 + 1.25331 I ``` --- Get results at high precision: ```wl In[1]:= NIntegrate[(1/1 + Sqrt[x]), {x, 0, 1}, PrecisionGoal -> 100, WorkingPrecision -> 200] Out[1]= 0.61370563888010938116553575708364686384899973127948949175863998101321275606061056878827334600716262491599703795885862853262895952848373888593465849672984807613854988384502718906646349188763555174624320 ``` --- Integrate along a complex line, finite or infinite: ```wl In[1]:= NIntegrate[Sqrt[x], {x, I, 3 - I}] Out[1]= 3.79214 - 2.21131 I In[2]:= NIntegrate[(Sinh[x]/x), {x, 0, I ∞}] Out[2]= -9.483638083064695`*^-14 + 1.5708 I ``` --- Along a piecewise linear contour in the complex plane: ```wl In[1]:= NIntegrate[(1/z + 1 / 2), {z, 1, E^(I π/3), E^(2 I π/3), -1, E^-(2 I π/3), E^-(I π/3), 1}] Out[1]= -1.1102230246251565`*^-16 + 6.28319 I ``` Along a circular contour in the complex plane: ```wl In[2]:= NIntegrate[(I Exp[I ω]/Exp[I ω] + 1 / 2), {ω, 0, 2π}] Out[2]= -1.1570744362643381`*^-12 + 6.28319 I ``` Plot the real part of a function and paths of integration: ```wl In[3]:= lpath = Graphics@Arrow@({Re[#], Im[#]}& /@ Exp[2π I / 6 Range[0, 6]]); In[4]:= cpath = ParametricPlot[{Cos[θ], Sin[θ]}, {θ, 0, 2π}, Axes -> None] /. Line -> Arrow; In[5]:= Show[DensityPlot[Re[(1/a + 1 / 2 + I b)], {a, -1.2, 1.2}, {b, -1.2, 1.2}], lpath, cpath] Out[5]= [image] ``` --- Vector- and tensor-valued functions: ```wl In[1]:= NIntegrate[{x, (1/Sqrt[x]), Sin[x]}, {x, 0, 5}] Out[1]= {12.5, 4.47214, 0.716338} In[2]:= NIntegrate[(⁠| | | | | ------- | ------- | ------- | | x | Sqrt[x] | x^1 / 3 | | Sqrt[x] | x^1 / 3 | x^1 / 4 | | x^1 / 3 | x^1 / 4 | x^1 / 5 |⁠), {x, 0, 5}]//MatrixForm Out[2]//MatrixForm= (⁠| | | | | ------- | ------- | ------- | | 12.5 | 7.45356 | 6.41241 | | 7.45356 | 6.41241 | 5.9814 | | 6.41241 | 5.9814 | 5.74887 |⁠) ``` --- Multivariate functions: ```wl In[1]:= NIntegrate[(BesselJ[3, y]/x + 1), {x, 0, 5}, {y, 0, 5}] Out[1]= 1.9431 In[2]:= NIntegrate[Exp[-x + y], {x, 0, ∞}, {y, 0, 1}] Out[2]= 1.71828 ``` --- The bounds for an integration variable can depend on the previous variables: ```wl In[1]:= NIntegrate[Sin[5x y + y ^ 2] + 1, {x, -1, 1}, {y, -Sqrt[1 - x ^ 2], Sqrt[1 - x ^ 2]}] Out[1]= 3.43678 ``` Plot the integrand over the integration region: ```wl In[2]:= Plot3D[Sin[5x y + y ^ 2] + 1, {x, -1, 1}, {y, -Sqrt[1 - x ^ 2], Sqrt[1 - x ^ 2]}, Filling -> 0, Mesh -> None] Out[2]= [image] ``` --- Integrate over regions defined by inequalities: ```wl In[1]:= NIntegrate[Boole[x ^ 2 + 2y ^ 2 + 3z ^ 2 < 1∧-1 < x < y < z < 1], {x, -∞, ∞}, {y, -∞, ∞}, {z, -∞, ∞}] Out[1]= 0.301328 ``` Plot the integration region: ```wl In[2]:= RegionPlot3D[x ^ 2 + 2y ^ 2 + 3z ^ 2 < 1∧-1 < x < y < z < 1, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None, PlotPoints -> 50] Out[2]= [image] ``` --- Integrate over any geometric region including special regions: ```wl In[1]:= NIntegrate[1, {x, y, z}∈Cone[]] Out[1]= 2.0944 In[2]:= Graphics3D[Cone[]] Out[2]= [image] ``` Formula regions: ```wl In[3]:= ℛ = ImplicitRegion[x ^ 2 + 2y ^ 2 + 3z ^ 2 < 1∧-1 < x < y < z < 1, {x, y, z}]; In[4]:= NIntegrate[1, {x, y, z}∈ℛ] Out[4]= 0.301328 ``` Mesh regions: ```wl In[5]:= ℛ = DelaunayMesh[RandomReal[1, {50, 2}]] Out[5]= [image] In[6]:= NIntegrate[1, {x, y}∈ℛ] Out[6]= 0.641144 ``` --- Find the volume of the unit ball in five dimensions: ```wl In[1]:= NIntegrate[1, {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3], Subscript[x, 4], Subscript[x, 5]}∈Ball[5]] Out[1]= 5.26379 ``` The surface area of the standard cylinder, by computing a surface integral: ```wl In[2]:= NIntegrate[1, {x, y, z}∈RegionBoundary[Cylinder[]]] Out[2]= 18.8496 ``` The length of an ellipse, by computing a line integral: ```wl In[3]:= NIntegrate[1, {x, y}∈Circle[{0, 0}, {2, 1}]] Out[3]= 9.68845 ``` --- Multivariate integrals using Monte Carlo methods: ```wl In[1]:= NIntegrate[1 / (1 + Underoverscript[∑, i, 100]x[i] ^ i), Evaluate[Sequence@@Table[{x[i], 0, 1}, {i, 100}]]] Out[1]= 0.224184 ``` #### Singular Integrals (9) Integrate functions with algebraic singularities at the endpoints: ```wl In[1]:= NIntegrate[{(1/Sqrt[x]), (1/x^1 / 5), (1/x^4 / 5)}, {x, 0, 1}] Out[1]= {2., 1.25, 5.} In[2]:= NIntegrate[{(1/Sqrt[1 - x]), (1/(1 - x)^1 / 5), (1/(1 - x)^4 / 5)}, {x, 0, 1}] Out[2]= {2., 1.25, 5.} ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, 0, 1}, PlotLabel -> Row[f, ", "], Filling -> Axis, PlotRange -> {0, 4}], {f, {{(1/Sqrt[x]), (1/x^1 / 5), (1/x^4 / 5)}, {(1/Sqrt[1 - x]), (1/(1 - x)^1 / 5), (1/(1 - x)^4 / 5)}}}] Out[3]= {[image], [image]} ``` --- Logarithmic singularities at endpoints: ```wl In[1]:= NIntegrate[{Log[x], Log[Log[x]]}, {x, 0, 1}] Out[1]= {-1., -0.577216 + 3.14159 I} In[2]:= NIntegrate[{Sqrt[Log[x]], Log[x] ^ 2}, {x, 0, 1}] Out[2]= {0. + 0.886227 I, 2.} ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, 0, 1}, Filling -> Axis, PlotLabel -> Row[f, ", "]], {f, {{Log[x], Re@Log[Log[x]]}, {Im@Sqrt[Log[x]], Log[x] ^ 2}}}] Out[3]= {[image], [image]} ``` --- Singularities at both endpoints: ```wl In[1]:= NIntegrate[{(1/Sqrt[1 - x] x^1 / 5), (1/Sqrt[1 - x] x^4 / 5), (Log[2 x]/Sqrt[1 - x])}, {x, 0, 1}] Out[1]= {2.29929, 6.26865, 0.158883} ``` Plot over the integration range: ```wl In[2]:= Plot[{(1/Sqrt[1 - x] x^1 / 5), (1/Sqrt[1 - x] x^4 / 5), (Log[2 x]/Sqrt[1 - x])}, {x, 0, 1}] Out[2]= [image] ``` --- Test for singularities interior to the integration region using ``Exclusions`` : ```wl In[1]:= NIntegrate[(1/Sqrt[Sin[x]]), {x, 0, 10}, Exclusions -> Sin[x] == 0] Out[1]= 10.4882 - 6.76947 I In[2]:= NIntegrate[Sqrt[Log[x - 7]], {x, 0, 20}, Exclusions -> {7, 8}] Out[2]= 25.8817 + 8.56065 I ``` Plot over the integration range: ```wl In[3]:= Table[Plot[{Re[f], Im[f]}, {x, 0, 20}, PlotLabel -> Style[{Re[f], Im[f]}, Small], Filling -> Axis], {f, {(1/Sqrt[Sin[x]]), Sqrt[Log[x - 7]]}}] Out[3]= [image] ``` --- Multivariate integrals of functions with singularities at edges and corners: ```wl In[1]:= NIntegrate[{(1/Sqrt[x]), Log[y], (Log[y]/x^4 / 5)}, {x, 0, 1}, {y, 0, 1}] Out[1]= {2., -1., -5.} In[2]:= NIntegrate[{(1/Sqrt[x + y]), Log[x y]}, {x, 0, 1}, {y, 0, 1}] Out[2]= {1.10457, -2.} ``` Plot over the integration range: ```wl In[3]:= Table[Plot3D[f, {x, 0, 1}, {y, 0, 1}, PlotLabel -> f, PlotPoints -> 35, Mesh -> None], {f, {(1/Sqrt[x]), Log[y], (Log[y]/x^4 / 5)}}] Out[3]= [image] ``` --- Singularities at the corners of multidimensional integration regions: ```wl In[1]:= NIntegrate[{(1/Sqrt[x + y]), Log[x + y], (Log[x + y]/Sqrt[1 - x y])}, {x, 0, 1}, {y, 0, 1}] Out[1]= {1.10457, -0.113706, -0.0341229} ``` Plot over the integration range: ```wl In[2]:= Table[Plot3D[f, {x, 0, 1}, {y, 0, 1}, PlotLabel -> f, Mesh -> None], {f, {(1/Sqrt[x + y]), Log[x + y], (Log[Sqrt[x + y]]/Sqrt[1 - x y])}}] Out[2]= {[image], [image], [image]} ``` --- Multivariate integral with a singularity at an interior point: ```wl In[1]:= NIntegrate[Log[2 - Sin[x] - Sin[y]], {x, 0, 4}, {y, 0, 4}, Exclusions -> {{(π/2), (π/2)}}] Out[1]= -2.41908 ``` Several interior point singularities: ```wl In[2]:= pl = RandomReal[{0, 4}, {10, 2}]; f = Sum[1 / Norm[{x, y} - p], {p, pl}]; In[3]:= NIntegrate[f, {x, 0, 4}, {y, 0, 4}, Exclusions -> pl]//Quiet Out[3]= 120.57 ``` Plot over the integration range: ```wl In[4]:= Table[Plot3D[g, {x, 0, 4}, {y, 0, 4}, Mesh -> None, PlotStyle -> Yellow], {g, {Log[2 - Sin[x] - Sin[y]], f}}] Out[4]= {[image], [image]} ``` --- Multivariate integral of a function that is singular at curves satisfying an equation: ```wl In[1]:= NIntegrate[(1/Sqrt[Cos[x y]]), {x, 0, 5}, {y, 0, 5}, Exclusions -> Cos[x y] == 0] Out[1]= 21.5381  - 19.5357 I ``` Plot over the integration range: ```wl In[2]:= Plot3D[{Re[(1/Sqrt[Cos[x y]])], Im[(1/Sqrt[Cos[x y]])]}, {x, 0, 5}, {y, 0, 5}, Exclusions -> Cos[x y] == 0, Mesh -> None, PlotStyle -> {Yellow, Magenta}, PlotPoints -> 35] Out[2]= [image] ``` --- Cauchy principal value of functions with non-integrable singularities: ```wl In[1]:= NIntegrate[(1/x - x^2), {x, -1, 2}, Method -> "PrincipalValue", Exclusions -> x - x^2 == 0] Out[1]= 1.38629 In[2]:= NIntegrate[(1/Log[x]), {x, 0, 10}, Method -> "PrincipalValue", Exclusions -> 1] Out[2]= 6.1656 ``` #### Piecewise Integrals (9) Piecewise functions with finitely many cases: ```wl In[1]:= NIntegrate[UnitStep[(x - 1)(x - 2)(x - 5)(x - 6)], {x, 0, 10}] Out[1]= 8. In[2]:= NIntegrate[Boole[(x - 1)(x - 2) ≤ (x - 5)(x - 6)], {x, 0, 10}] Out[2]= 3.5 In[3]:= NIntegrate[Min[(x - 1)(x - 3), (x - 5)(x - 7)], {x, 0, 10}] Out[3]= 19.3333 In[4]:= Table[Plot[f, {x, 0, 10}, Filling -> Axis, PlotLabel -> Style[f, Small]], {f, {UnitStep[(x - 1)(x - 2)(x - 5)(x - 6)], Boole[(x - 1)(x - 2) ≤ (x - 5)(x - 6)], Min[(x - 1)(x - 3), (x - 5)(x - 7)]}}] Out[4]= {[image], [image], [image]} ``` --- Piecewise functions with infinitely many cases: ```wl In[1]:= NIntegrate[Floor[x ^ 2], {x, 0, 5}] Out[1]= 39.3662 In[2]:= NIntegrate[FractionalPart[x], {x, 0, 5}] Out[2]= 2.5 In[3]:= NIntegrate[PrimePi[x ^ 2], {x, 0, 5}] Out[3]= 16.7719 In[4]:= Table[Plot[f, {x, 0, 5}, Filling -> Axis, PlotLabel -> Style[f, Small], PlotPoints -> 100], {f, {Floor[x ^ 2], FractionalPart[x], PrimePi[x ^ 2]}}] Out[4]= {[image], [image], [image]} ``` --- Composition with any other function: ```wl In[1]:= NIntegrate[Exp[UnitStep[(x - 1)(x - 2)(x - 5)]], {x, 0, 10}] Out[1]= 20.3097 In[2]:= NIntegrate[Sin[Min[(x - 1), (x - 5)(x - 7)]], {x, 0, 10}] Out[2]= 2.16896 In[3]:= NIntegrate[BesselJ[2, FractionalPart[x]], {x, 0, 10}] Out[3]= 0.396292 In[4]:= Table[Plot[f, {x, 0, 5}, Filling -> Axis, PlotLabel -> Style[f, Small], PlotPoints -> 100], {f, {Exp[UnitStep[(x - 1)(x - 2)(x - 5)]], Sin[Min[(x - 1), (x - 5)(x - 7)]], BesselJ[2, FractionalPart[x]]}}] Out[4]= {[image], [image], [image]} ``` --- Use ``Exclusions`` to explicitly specify discontinuities or sharp corners: ```wl In[1]:= NIntegrate[{Log[-3 + I x], Sqrt[-1 + I x]}, {x, -2, 3}, Exclusions -> {0}] Out[1]= {6.02071 + 2.44954 I, 2.65134 + 1.35829 I} In[2]:= NIntegrate[{Re[EllipticF[x, 2]], Im[EllipticF[x, 2]]}, {x, -2, 3}, Exclusions -> Csc[x]^2 == 2] Out[2]= {1.8041, -2.42089} ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, -2, 3}, Filling -> Axis, PlotLabel -> Row[f, ", "]], {f, {{Im[Log[-3 + I x]], Im[Sqrt[-1 + I x]]}, {Re[EllipticF[x, 2]], Im[EllipticF[x, 2]]}}}] Out[3]= {[image], [image]} ``` --- Multivariate piecewise integrals with finitely many cases: ```wl In[1]:= NIntegrate[Max[ x y^2, x^2y], {x, -2, 2}, {y, -2, 2}] Out[1]= 12.8 In[2]:= NIntegrate[Max[Abs[x ^ 2 - y], Abs[x - y ^ 2]], {x, -2, 2}, {y, -2, 2}] Out[2]= 35.7677 - 3.269229368127778`*^-32 I In[3]:= NIntegrate[x y Floor[x y], {x, -2, 2}, {y, -2, 2}] Out[3]= 28.504 In[4]:= Table[Plot3D[f, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 35, Mesh -> None, ExclusionsStyle -> {None, Red}, PlotLabel -> f], {f, {Max[ x y^2, x^2y], Max[ Abs[ x ^ 2 - y], Abs[x - y ^ 2]], x y Floor[x y]}}] Out[4]= [image] ``` --- Explicitly specify discontinuities and sharp corners: ```wl In[1]:= NIntegrate[Sqrt[x + I y], {x, -2, 2}, {y, -2, 2}, Exclusions -> y == 0] Out[1]= 12.3924 - 2.6645352591003757`*^-15 I ``` Plot real and imaginary parts over the region: ```wl In[2]:= Plot3D[{Re[Sqrt[x + I y]], Im[Sqrt[x + I y]]}, {x, -2, 2}, {y, -2, 2}, Mesh -> None, ExclusionsStyle -> {None, Red}, PlotStyle -> {Lighter[Purple, 0.5], White}] Out[2]= [image] ``` --- Integrate over 2D regions: ```wl In[1]:= NIntegrate[Boole[x ≤ 2y + 1 && y > x - 1], {x, -2, 2}, {y, -2, 2}] Out[1]= 9.75 In[2]:= NIntegrate[Boole[x ^ 2 + y ^ 2 < 1], {x, -2, 2}, {y, -2, 2}] Out[2]= 3.14159 In[3]:= NIntegrate[Boole[1 / 4 ≤ x ^ 2 + y ^ 2 ≤ 1 && y ≥ 0 && y ≥ -x], {x, -2, 2}, {y, -2, 2}] Out[3]= 0.883573 In[4]:= Table[RegionPlot[r, {x, -2, 2}, {y, -2, 2}, PlotLabel -> Style[r, Small]], {r, {x ≤ 2y + 1 && y > x - 1, x ^ 2 + y ^ 2 < 1, 1 / 4 ≤ x ^ 2 + y ^ 2 ≤ 1 && y ≥ 0 && y ≥ -x}}] Out[4]= {[image], [image], [image]} ``` --- Integrate over 3D regions: ```wl In[1]:= NIntegrate[Boole[x ≤ 2y + 3z + 1 && y > x - z], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[1]= 29.8611 In[2]:= NIntegrate[Boole[x ^ 2 + 2y ^ 2 + z ^ 2 ≤ 1], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[2]= 2.96192 In[3]:= NIntegrate[(2x^2 + 3y^2) Boole[x^2 + y^2 + z^2 ≤ 1∧3x^2 + 3y^2 ≤ z^2], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}] Out[3]= 0.107742 In[4]:= Table[RegionPlot3D[r, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, PlotLabel -> r], {r, {x ≤ 2y + 3z + 1 && y > x - z, x ^ 2 + 2y ^ 2 + z ^ 2 ≤ 1, x^2 + y^2 + z^2 ≤ 1∧3x^2 + 3y^2 ≤ z^2}}] Out[4]= [image] ``` --- Integrate over $n$-dimensional regions: ```wl In[1]:= NIntegrate[Boole[x ≤ 2y + 3z + w + 1 && y > x - z + w], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, {w, -2, 2}] Out[1]= 101.42 In[2]:= NIntegrate[Boole[x ^ 2 + y ^ 2 + z ^ 2 + w ^ 2 ≤ 1], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, {w, -2, 2}] Out[2]= 4.9348 ``` #### Oscillatory Integrals (15) Integrating a highly oscillatory elementary function over a finite range: ```wl In[1]:= NIntegrate[{Cos[x], Sin[x]}, {x, 1, 10 ^ 4}] Out[1]= {-1.14709, 1.49246} In[2]:= NIntegrate[{Exp[x I], 2 ^ (x I)}, {x, 1, 10 ^ 4}] Out[2]= {-1.14709 + 1.49246 I, 0.375744 + 0.479159 I} ``` Plot over $\frac{1}{100}$ the range: ```wl In[3]:= Table[Plot[f, {x, 1, 100}, PlotLabel -> f], {f, {{Cos[x], Sin[x]}, Im@{Exp[x I], 2 ^ (x I)}}}] Out[3]= {[image], [image]} ``` --- Highly oscillatory special functions: ```wl In[1]:= NIntegrate[{BesselJ[1, x], BesselY[2, x]}, {x, 1, 10 ^ 4}] Out[1]= {0.772294, -0.932453} In[2]:= NIntegrate[{BesselI[1, x I], BesselK[2, x I]}, {x, 1, 10 ^ 4}] Out[2]= {1.469487575124191`*^-17 + 0.772294 I, -1.46469 + 1.50282 I} In[3]:= NIntegrate[{AiryAi[-x], AiryBi[-x] / Sqrt[x]}, {x, 0, 10 ^ 4}] Out[3]= {0.667162, 0.628585} In[4]:= NIntegrate[{SinIntegral[x], CosIntegral[x]}, {x, 1, 10 ^ 4}] Out[4]= {15706.5, 0.504162} In[5]:= NIntegrate[{FresnelC[x / 10], FresnelS[x / 10]}, {x, 0, 10 ^ 4}] Out[5]= {5000., 4996.82} ``` Plot over $\frac{1}{100}$ the range: ```wl In[6]:= Table[Plot[f, {x, 1, 100}, PlotLabel -> Row[f, ","]], {f, {{BesselJ[1, x], BesselY[2, x]}, Im@{BesselI[1, x I], BesselK[2, x I]}, {AiryAi[-x], AiryBi[-x] / Sqrt[x]}, {SinIntegral[x], CosIntegral[x]}, {FresnelC[x / 10], FresnelS[x / 10]}}}] Out[6]= [image] ``` --- Integrate oscillatory functions over an infinite range: ```wl In[1]:= NIntegrate[{Cos[x], Sin[x]}Exp[-x ^ 2], {x, 1, ∞}] Out[1]= {0.0340199, 0.129738} In[2]:= NIntegrate[{Exp[x I], 2 ^ (x I)}Exp[-x ^ 2], {x, 1, ∞}] Out[2]= {0.0340199 + 0.129738 I, 0.0836617 + 0.108261 I} ``` Oscillatory special functions over an infinite range: ```wl In[3]:= NIntegrate[{BesselJ[1, x], BesselY[2, x]}, {x, 1, ∞}] Out[3]= {0.765198, -0.925356} In[4]:= NIntegrate[{BesselI[1, x I], BesselK[2, x I]}, {x, 1, ∞}] Out[4]= {0. + 0.765198 I, -1.45355 + 1.50855 I} In[5]:= NIntegrate[{AiryAi[-x], AiryBi[-x]}, {x, 1, ∞}] Out[5]= {0.200993, -0.373005} ``` --- Sums of oscillatory functions: ```wl In[1]:= NIntegrate[Sin[Sqrt[2]x] + Sin[Sqrt[3]x] + Sin[Sqrt[5]x], {x, 1, 10 ^ 4}] Out[1]= -0.254076 In[2]:= NIntegrate[BesselJ[1, x] + FresnelS[x / 10], {x, 1, 10 ^ 4}] Out[2]= 4997.59 ``` Plot over $\frac{1}{100}$ the range: ```wl In[3]:= Table[Plot[f, {x, 1, 100}, PlotLabel -> Style[f, Small]], {f, {Sin[Sqrt[2]x] + Sin[Sqrt[3]x] + Sin[Sqrt[5]x], BesselJ[1, x] + FresnelS[x / 10]}}] Out[3]= {[image], [image]} ``` --- Products of oscillatory functions: ```wl In[1]:= NIntegrate[Sin[Sqrt[2]x]Sin[Sqrt[3]x]Sin[Sqrt[5]x], {x, 1, 10 ^ 4}] Out[1]= 0.134673 In[2]:= NIntegrate[BesselJ[1, x]BesselJ[2, Sqrt[2]x], {x, 1, 10 ^ 4}] Out[2]= 0.474294 ``` Plot over $\frac{1}{100}$ the range: ```wl In[3]:= Table[Plot[f, {x, 1, 100}, PlotLabel -> Style[f, Small]], {f, {Sin[Sqrt[2]x]Sin[Sqrt[3]x]Sin[Sqrt[5]x], BesselJ[1, x]BesselJ[2, Sqrt[2]x]}}] Out[3]= {[image], [image]} ``` --- Powers of oscillatory functions: ```wl In[1]:= NIntegrate[BesselJ[2, Sqrt[2]x]^2, {x, 1, 10 ^ 4}] Out[1]= 2.13898 In[2]:= NIntegrate[(Sin[Sqrt[2]x]Sin[Sqrt[3]x])^3, {x, 1, 10 ^ 4}] Out[2]= -1.04541 ``` Plot over $\frac{1}{100}$ the range: ```wl In[3]:= Table[Plot[f, {x, 1, 100}, PlotLabel -> Style[f, Small]], {f, {BesselJ[2, Sqrt[2]x]^2, (Sin[Sqrt[2]x]Sin[Sqrt[3]x])^3}}] Out[3]= {[image], [image]} ``` --- Compositions of oscillatory functions with non-oscillatory ones: ```wl In[1]:= NIntegrate[Sin[(1/(5x - 1)(5x - 2)(5x - 3))], {x, 0, 1}] Out[1]= 0.0169585 In[2]:= NIntegrate[BesselJ[3 / 7, (Log[x]/x)], {x, 0, 1}] Out[2]= 0.0660302 + 0.289297 I ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, 0, 1}, PlotLabel -> f], {f, {Sin[(1/(5x - 1)(5x - 2)(5x - 3))], Re[BesselJ[3 / 7, (Log[x]/x)]]}}] Out[3]= {[image], [image]} ``` --- Sums, products, powers, and compositions of oscillatory functions: ```wl In[1]:= NIntegrate[BesselY[2 / 3, x]Sin[x + Sqrt[x]] ^ 4, {x, 0, 10 ^ 4}] Out[1]= -0.387051 In[2]:= NIntegrate[Sin[x] (Sin[x ^ 2] (Sin[x ^ 3] + (1/x)) + (1/x)), {x, 0, ∞}] Out[2]= 2.43879 ``` Plot over part of the integration range: ```wl In[3]:= Table[Plot[f, {x, 0, 20}, PlotLabel -> Style[f, Small]], {f, {BesselY[2 / 3, x]Sin[x + Sqrt[x]] ^ 4, Sin[x] ((1/x) + Sin[x ^ 2] ((1/x) + Sin[x ^ 3]))}}] Out[3]= {[image], [image]} ``` --- Multivariate integrals of highly oscillatory elementary functions over a finite range: ```wl In[1]:= NIntegrate[{Sin[x]Cos[y], Sin[x + y]}, {x, 1, 10 ^ 4}, {y, 1, 10 ^ 4}] Out[1]= {-1.71198, -3.42395} In[2]:= NIntegrate[{Exp[I(x + y)], 2 ^ (I x y)}, {x, 1, 10 ^ 4}, {y, 1, 10 ^ 4}] Out[2]= {-0.911625 - 3.42395 I, -1.29231 - 0.133751 I} ``` Plot over one millionth of the range: ```wl In[3]:= ParallelTable[Plot3D[f, {x, 1, 10}, {y, 1, 10}, PlotLabel -> f, Mesh -> False, PlotPoints -> 50], {f, {Sin[x]Cos[y], Sin[x + y], Re[Exp[I(x + y)]], Re[2 ^ (I x y)]}}] Out[3]= [image] ``` --- Multivariate integrals of oscillatory special functions over a finite range: ```wl In[1]:= NIntegrate[{BesselJ[1, x]BesselY[2, y], BesselJ[3, x + y]}, {x, 10 ^ 4, 2 10 ^ 4}, {y, 10 ^ 4, 2 10 ^ 4}] Out[1]= {-0.000160337, 0.00218355} In[2]:= NIntegrate[{BesselI[1, x I]BesselK[2, y I], BesselI[3, I(x + y)]}, {x, 10 ^ 4, 2 10 ^ 4}, {y, 10 ^ 4, 2 10 ^ 4}] Out[2]= {0.0000908712 - 0.000251858 I, -5.038725795352576`*^-17 - 0.00218355 I} ``` --- Multivariate oscillatory functions over an infinite range: ```wl In[1]:= NIntegrate[{(Sin[x + y]/1 + x ^ 4 + y ^ 4), (Sin[x] + Cos[y]/1 + x ^ 4 + y ^ 4)}, {x, 0, ∞}, {y, 0, ∞}] Out[1]= {0.88768, 1.60749} In[2]:= NIntegrate[Exp[I(x ^ 2 + y ^ 2)], {x, -∞, ∞}, {y, -∞, ∞}, PrecisionGoal -> 5]//Quiet Out[2]= -1.3458767988971942`*^-9 + 3.14159 I In[3]:= NIntegrate[(BesselY[1, x + y]/x ^ 4 + y ^ 4), {x, 1, ∞}, {y, 1, ∞}] Out[3]= 0.042943 In[4]:= NIntegrate[Sin[Exp[x] + y ^ 2], {x, 0, ∞}, {y, 0, ∞}] Out[4]= 0.47094 ``` --- Singular oscillatory functions: ```wl In[1]:= NIntegrate[{(1/Sqrt[x]) + Cos[x], (1/Sqrt[x]) + BesselY[0, x]}, {x, 0, 5000}] Out[1]= {140.433, 141.428} In[2]:= NIntegrate[{(Sin[x]/Sqrt[x]), (BesselJ[1, x]/Sqrt[x])}, {x, 0, 5000}] Out[2]= {1.25113, 0.956072} ``` Plot over $\frac{1}{10}$ the range: ```wl In[3]:= Table[Plot[f, {x, 0, 500}, PlotLabel -> f], {f, {{(1/Sqrt[x]) + Cos[x], (1/Sqrt[x]) + BesselY[0, x]}, {(Cos[x]/Sqrt[x]), (BesselJ[1, x]/Sqrt[x])}}}] Out[3]= [image] ``` --- Composition of oscillatory functions with singular ones: ```wl In[1]:= NIntegrate[{Sin[10 ^ 3Sqrt[x]], AiryAi[10Log[x]]}, {x, 0, 1}] Out[1]= {-0.0011231, 0.0639225} In[2]:= NIntegrate[{Sin[(1/x)]Cos[(1/1 - x)], (BesselJ[1, 1 / x]/x)}, {x, 0, 1}] Out[2]= {-0.210282, 0.52032} ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, 0, 1}, PlotLabel -> f], {f, {{Sin[10 ^ 3Sqrt[x]], AiryAi[10Log[x]]}, {Sin[(1/x)]Cos[(1/1 - x)], (BesselJ[1, 1 / x]/x)}}}] Out[3]= {[image], [image]} ``` --- Piecewise oscillatory functions: ```wl In[1]:= NIntegrate[FractionalPart[x ^ 2]Sin[10 ^ 4x], {x, 0, 3}] Out[1]= 0.000202256 In[2]:= NIntegrate[Round[x] AiryAi[-10 ^ 4 TriangleWave[x] ^ 2], {x, 0, 3}] Out[2]= 0.0244967 ``` Plot at $\frac{1}{100}$ the oscillation rate: ```wl In[3]:= Table[Plot[f, {x, 0, 3}, PlotLabel -> f], {f, {FractionalPart[x ^ 2]Sin[10 ^ 2x], x AiryAi[-10 ^ 2TriangleWave[x] ^ 2]}}] Out[3]= {[image], [image]} ``` --- Piecewise oscillatory functions with singularities: ```wl In[1]:= NIntegrate[Cos[(x/FractionalPart[x])], {x, -2, 3}] Out[1]= -0.0996593 In[2]:= NIntegrate[(Sin[10 ^ 2x ^ 4]/Sqrt[x - Floor[x]]), {x, -2, 3}, MaxRecursion -> 20] Out[2]= 0.237224 ``` Plot over the integration range: ```wl In[3]:= Table[Plot[f, {x, -2, 3}, PlotLabel -> f], {f, {Cos[(x/FractionalPart[x])], (Sin[10 ^ 2x ^ 4]/Sqrt[x - Floor[x]])}}] Out[3]= {[image], [image]} ``` #### Region Integrals (11) Different ways of integrating over a unit disk: ```wl In[1]:= NIntegrate[Boole[x ^ 2 + y ^ 2 < 1], {x, -1, 1}, {y, -1, 1}] Out[1]= 3.14159 In[2]:= NIntegrate[1, {x, y}∈Disk[]] Out[2]= 3.14159 ``` Visualize the integrand over the region: ```wl In[3]:= Plot3D[{1, 0}, {x, y}∈Disk[]] Out[3]= [image] ``` --- More general integral over a unit disk: ```wl In[1]:= NIntegrate[Abs[Sqrt[x]]Boole[x ^ 2 + y ^ 2 < 1], {x, -1, 1}, {y, -1, 1}] Out[1]= 1.91702 In[2]:= NIntegrate[Abs[Sqrt[x]], {x, y}∈Disk[]] Out[2]= 1.91702 ``` Visualize the integrand over the region: ```wl In[3]:= Plot3D[{Abs[Sqrt[x]], 0}, {x, y}∈Disk[]] Out[3]= [image] ``` --- Integration over a 0D region corresponds to summing over the points: ```wl In[1]:= NIntegrate[x, x∈Point[{{2}, {3}}]] Out[1]= 5. In[2]:= NIntegrate[{x, y}, {x, y}∈Point[{{1, 2}, {3, 4}}]] Out[2]= {4., 6.} In[3]:= NIntegrate[{x, y, z}, {x, y, z}∈Point[{{1, 2, 3}, {4, 5, 6}}]] Out[3]= {5., 7., 9.} ``` --- Integration over a 1D region corresponds to a curve integral: ```wl In[1]:= NIntegrate[1, {x, y}∈Circle[]] Out[1]= 6.28319 In[2]:= NIntegrate[1, {x, y, z}∈ImplicitRegion[x ^ 2 + y ^ 2 == 1∧z == 0, {x, y, z}]] Out[2]= 6.28319 ``` --- Integration over a 2D region corresponds to a surface integral: ```wl In[1]:= NIntegrate[1, {x, y, z}∈Sphere[]] Out[1]= 12.5664 ``` --- In general, integrals over lower-dimensional regions correspond to surface integrals: ```wl In[1]:= NIntegrate[1, {x, y, z, w}∈Sphere[{0, 0, 0, 0}, 1]] Out[1]= 19.7392 ``` The region is lower dimensional: ```wl In[2]:= RegionDimension[Sphere[{0, 0, 0, 0}, 1]] < RegionEmbeddingDimension[Sphere[{0, 0, 0, 0}, 1]] Out[2]= True ``` --- Integrate over any basic region including a filled parallelogram: ```wl In[1]:= NIntegrate[Sin[x]Cos[y], {x, y}∈Parallelogram[]] Out[1]= 0.632261 In[2]:= Graphics[{LightBlue, Parallelogram[]}, ImageSize -> Tiny] Out[2]= [image] ``` A filled cone: ```wl In[3]:= NIntegrate[x ^2 + y z , {x, y, z}∈Cone[]] Out[3]= 0.314159 In[4]:= Graphics3D[Cone[], ImageSize -> Tiny] Out[4]= [image] ``` A sphere surface: ```wl In[5]:= NIntegrate[x ^2 + Sin[y z ], {x, y, z}∈Sphere[]] Out[5]= 4.18879 In[6]:= Graphics3D[Sphere[], ImageSize -> Tiny] Out[6]= [image] ``` A filled ellipsoid: ```wl In[7]:= NIntegrate[x ^2 + y^2 + z^2 , {x, y, z}∈Ellipsoid[{0, 0, 0}, {1, 2, 3}]] Out[7]= 70.3717 In[8]:= Graphics3D[Ellipsoid[{0, 0, 0}, {1, 2, 3}], ImageSize -> Tiny] Out[8]= [image] ``` A half-infinite line: ```wl In[9]:= NIntegrate[Exp[-x^2]y, {x, y}∈HalfLine[{0, 0}, {1, 1}]] Out[9]= 0.707107 In[10]:= Graphics[{HalfLine[{0, 0}, {1, 1}]}, ImageSize -> Tiny] Out[10]= [image] ``` --- Integrate over any ``ImplicitRegion``: ```wl In[1]:= ℛ = ImplicitRegion[x^2 + y^2 == 1, {x, y}]; In[2]:= NIntegrate[ x ^ 2 + y ^ 2, {x, y}∈ℛ] Out[2]= 6.28319 In[3]:= RegionPlot[ℛ, Frame -> False, ImageSize -> Tiny] Out[3]= [image] ``` A full-dimensional region: ```wl In[4]:= ℛ = ImplicitRegion[x ^ 3 + y ^ 3 - x y ≤ 1 && x ≥ 0 && y ≥ -2, {x, y}]; In[5]:= NIntegrate[1, {x, y}∈ℛ] Out[5]= 3.45155 In[6]:= RegionPlot[ℛ, Frame -> False, ImageSize -> Tiny] Out[6]= [image] ``` --- Integrate over any ``ParametricRegion``: ```wl In[1]:= ℛ = ParametricRegion[{Cos[θ], Sin[θ]}, {{θ, 0, 2π}}]; In[2]:= NIntegrate[x ^ 2 + y ^ 2, {x, y}∈ℛ] Out[2]= 6.28319 In[3]:= RegionPlot[ℛ, Frame -> False, ImageSize -> Tiny] Out[3]= [image] ``` A full-dimensional region: ```wl In[4]:= ℛ = ParametricRegion[{{s, t s ^ 2}, s ^ 2 + t ^ 2 ≤ 1}, {s, t}]; In[5]:= NIntegrate[x ^ 2 + y ^ 2, {x, y}∈ℛ] Out[5]= 0.417243 In[6]:= RegionPlot[ℛ, Frame -> False, ImageSize -> Tiny] Out[6]= [image] ``` --- Integrate over any ``MeshRegion``: ```wl In[1]:= NIntegrate[x, x∈[image]] Out[1]= 3. ``` Regions in 2D: ```wl In[2]:= NIntegrate[x, {x, y}∈[image]] Out[2]= 2.82843 In[3]:= NIntegrate[x, {x, y}∈[image]] Out[3]= 3. ``` Regions in 3D: ```wl In[4]:= NIntegrate[x ^ 2 + y z , {x, y, z}∈[image]] Out[4]= 0.727083 In[5]:= NIntegrate[x ^ 2 + y ^ 2 + z ^ 2 , {x, y, z}∈[image]] Out[5]= 26.6667 ``` --- Integrate over any ``BoundaryMeshRegion`` : ```wl In[1]:= NIntegrate[x, {x, y}∈[image]] Out[1]= 12. In[2]:= NIntegrate[x ^ 2 + y z , {x, y, z}∈[image]] Out[2]= 4.6 ``` ### Options (30) #### AccuracyGoal (1) The ``AccuracyGoal`` option can be used to change the default absolute tolerance: ```wl In[1]:= exact = Integrate[10 ^ (-5) / (1 + 10 ^ 2(x - 1 / 2) ^ 2), {x, 0, 1}] Out[1]= (ArcTan[5]/500000) ``` The integration process stops once the accuracy goal criterion has been exceeded: ```wl In[2]:= NIntegrate[10 ^ (-5) / (1 + 10 ^ 2(x - 1 / 2) ^ 2), {x, 0, 1}, AccuracyGoal -> 8] - exact Out[2]= -4.9579365563708815`*^-12 ``` The result with the default settings is different since the default uses only a precision criterion: ```wl In[3]:= NIntegrate[10 ^ (-5) / (1 + 10 ^ 2(x - 1 / 2) ^ 2), {x, 0, 1}] - exact Out[3]= 4.466404730871926`*^-18 ``` #### EvaluationMonitor (2) Get the number of evaluation points used in a numerical integration: ```wl In[1]:= Block[{k = 0}, {NIntegrate[(1/Sqrt[x]), {x, 0, 1}, EvaluationMonitor :> k++], k}] Out[1]= {2., 132} ``` --- Show the evaluation points used in a numerical integration: ```wl In[1]:= ListPlot[Reap[NIntegrate[1 / Sqrt[x], {x, -1, 0, 1}, EvaluationMonitor :> Sow[x]]][[2, 1]], PlotRange -> All] Out[1]= [image] ``` #### Exclusions (1) Integration by excluding the curves on which the integrand's denominator is zero: ```wl In[1]:= NIntegrate[(1/Sqrt[Sin[x^2 + y]]), {x, -5, 5}, {y, -5, 5}, Exclusions -> (Sin[x^2 + y] == 0)] Out[1]= 81.8469  - 84.7051 I ``` The curves on which the integrand is singular: ```wl In[2]:= ContourPlot[Sin[x ^ 2 + y] == 0, {x, -5, 5}, {y, -5, 5}] Out[2]= [image] ``` #### MaxPoints (1) Stop integration after a specified number of points has been exceeded: ```wl In[1]:= NIntegrate[1 / Sqrt[x], {x, 0, 1}, MaxPoints -> 20] ``` NIntegrate::maxp: The integral failed to converge after 33 integrand evaluations. NIntegrate obtained 1.9558072180392028\` and 0.06781302015519788\` for the integral and error estimates. ```wl Out[1]= 1.95581 In[2]:= NIntegrate[1 / Sqrt[Abs[x + y + z]], {x, -1, 1}, {y, -1, 10}, {z, -1, 30}, Method -> "AdaptiveMonteCarlo", MaxPoints -> 1000] ``` NIntegrate::maxp: The integral failed to converge after 1100 integrand evaluations. NIntegrate obtained 197.45748669654552\` and 2.6568252158578005\` for the integral and error estimates. ```wl Out[2]= 197.457 ``` #### MaxRecursion (1) Without enough adaptive recursion, the following gives a poor result: ```wl In[1]:= NIntegrate[1 / Sqrt[Sin[x]], {x, 0, 10}] ``` NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {9.43625271977563}. NIntegrate obtained 10.3032-6.74255 I and 0.3168981771605228\` for the integral and error estimates. ```wl Out[1]= 10.3032 - 6.74255 I ``` Specifying a larger value for ``MaxRecursion`` gives a much better result: ```wl In[2]:= NIntegrate[1 / Sqrt[Sin[x]], {x, 0, 10}, MaxRecursion -> 100] ``` NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. ```wl Out[2]= 10.4882 - 6.76947 I ``` Specifying the singularity locations is even more efficient: ```wl In[3]:= NIntegrate[1 / Sqrt[Sin[x]], {x, 0, π, 2 π, 3 π, 10}] Out[3]= 10.4882 - 6.76947 I ``` #### Method (21) ##### Integration Rules (10) Left- and right-sided Riemann sum: ```wl In[2]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"RiemannRule", "Type" -> "Left"}, PrecisionGoal -> 2] Out[2]= 12.1732 In[3]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"RiemannRule", "Type" -> "Right"}, PrecisionGoal -> 2] Out[3]= 12.0847 ``` Riemann sum samples uniformly at the left or right endpoints of subregions: ```wl In[2]:= Table[ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"RiemannRule", "Type" -> t, "Points" -> 5}, MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, PlotLabel -> t, PlotRange -> {{-1, 11}}, Frame -> True, Axes -> False], {t, {"Left", "Right"}}]//Quiet Out[2]= {[image], [image]} ``` --- Basic trapezoidal rule with no extrapolation, corresponding to piecewise linear approximation: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"TrapezoidalRule", "RombergQuadrature" -> False}] Out[1]= 12.1561 ``` Trapezoidal rule with Romberg extrapolation: ```wl In[2]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "TrapezoidalRule"] Out[2]= 12.1561 ``` The basic trapezoidal rule samples uniformly (when adaptivity is turned off): ```wl In[3]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"TrapezoidalRule", "RombergQuadrature" -> False}, MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0]//Quiet Out[3]= [image] ``` Default adaptive method with trapezoidal rule: ```wl In[4]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"TrapezoidalRule", "RombergQuadrature" -> False}, MaxRecursion -> 5, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0]//Quiet Out[4]= [image] ``` --- Newton–Cotes rule with evenly spaced sampling points: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "NewtonCotesRule"] Out[1]= 12.1561 ``` Closed formulas include endpoints, but open formulas do not: ```wl In[2]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"NewtonCotesRule", "Type" -> "Closed"}] Out[2]= 12.1561 In[3]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"NewtonCotesRule", "Type" -> "Open"}] Out[3]= 12.1561 ``` A Newton–Cotes rule corresponds to polynomial interpolation: ```wl In[4]:= data = Last@Last@Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "NewtonCotesRule", MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]//Quiet Out[4]= {{0., 2.71828}, {2.5, 0.448815}, {5., 1.32798}, {7.5, 1.4143}, {10., 0.432112}} In[5]:= poly = Expand@InterpolatingPolynomial[data, x] Out[5]= 2.71828  - 2.42963 x + 0.836038 x^2 - 0.100696 x^3 + 0.00391022 x^4 In[6]:= Show[ListPlot[data, Filling -> Axis], Plot[{poly, Exp[Cos[x]]}, {x, 0, 10}]] Out[6]= [image] ``` The approximation with no adaptivity is the same as integrating the corresponding polynomial: ```wl In[7]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "NewtonCotesRule", MaxRecursion -> 0] ``` NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 0 recursive bisections in x near {x} = {10.}. NIntegrate obtained 10.845364976464321\` and 3.258518921949677\` for the integral and error estimates. ```wl Out[7]= 10.8454 In[8]:= Integrate[poly, {x, 0, 10}] Out[8]= 10.8454 ``` The method with order $n$ is exact for polynomials up to degree $n$: ```wl In[9]:= NIntegrate[1 + x + x ^ 3, {x, 0, 10}, Method -> {"NewtonCotesRule", "Order" -> 3}, MaxRecursion -> 0]//Quiet Out[9]= 2560. In[10]:= Integrate[1 + x + x ^ 3, {x, 0, 10}] Out[10]= 2560 ``` --- Clenshaw–Curtis quadrature rule with the strategy selected automatically: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "ClenshawCurtisRule"] Out[1]= 12.1561 ``` The sampling points are nonuniform: ```wl In[2]:= data = Last@Last@Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "ClenshawCurtisRule", MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]//Quiet Out[2]= {{0., 2.71828}, {0.380602, 2.53056}, {1.46447, 1.11197}, {3.08658, 0.368436}, {5., 1.32798}, {6.91342, 2.24317}, {8.53553, 0.532592}, {9.6194, 0.374891}, {10., 0.432112}} In[3]:= ListPlot[data, Filling -> Axis] Out[3]= [image] ``` The points are rescaled versions of $\cos (\pi k/8)$: ```wl In[4]:= 5. - 5Cos[Pi / 8 Range[0, 8]] Out[4]= {0., 0.380602, 1.46447, 3.08658, 5., 6.91342, 8.53553, 9.6194, 10.} In[5]:= First /@ data Out[5]= {0., 0.380602, 1.46447, 3.08658, 5., 6.91342, 8.53553, 9.6194, 10.} ``` --- Gaussian quadrature rule with Kronrod extension for error estimation: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "GaussKronrodRule"] Out[1]= 12.1561 ``` Gaussian rules use nonuniform sample points: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "GaussKronrodRule", MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0]//Quiet Out[2]= [image] ``` The method with $n$ Gauss points is exact for polynomials of degree $3 n+1$: ```wl In[3]:= Table[NIntegrate[x ^ (3n + 1), {x, 0, 2}, Method -> {"GaussKronrodRule", "Points" -> n}, MaxRecursion -> 0]//Quiet, {n, 5, 10}] Out[3]= {7710.12, 52428.8, 364722., 2.581110153846171*^6, 1.851279006896558*^7, 1.3421772799999928*^8} In[4]:= N@Table[Integrate[x ^ (3n + 1), {x, 0, 2}], {n, 5, 10}] Out[4]= {7710.12, 52428.8, 364722., 2.581110153846154*^6, 1.8512790068965517*^7, 1.34217728*^8} ``` The method with order $p$ uses enough points to be exact for polynomials of degree $p$ : ```wl In[5]:= Table[NIntegrate[x ^ p, {x, 0, 2}, Method -> {"GaussKronrodRule", "Order" -> p}, MaxRecursion -> 0]//Quiet, {p, 10, 15}] Out[5]= {186.182, 341.333, 630.154, 1170.29, 2184.53, 4096.} In[6]:= N@Table[Integrate[x ^ p, {x, 0, 2}], {p, 10, 15}] Out[6]= {186.182, 341.333, 630.154, 1170.29, 2184.53, 4096.} ``` --- Gaussian quadrature rule at Lobatto points with Kronrod extension: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "LobattoKronrodRule"] Out[1]= 12.1561 ``` Lobatto points are nonuniform and include the endpoints of the integration region: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "LobattoKronrodRule", MaxRecursion -> 0, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0]//Quiet Out[2]= [image] ``` --- Multipanel rule (or composite rule) applies the specified rule to multiple subintervals: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"MultipanelRule", Method -> "GaussKronrodRule", "Panels" -> 8}] Out[1]= 12.1561 ``` Any other rule can be used together with a multipanel rule: ```wl In[2]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"MultipanelRule", Method -> "NewtonCotesRule", "Panels" -> 8}] Out[2]= 12.1561 In[3]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> {"MultipanelRule", Method -> "ClenshawCurtisRule", "Panels" -> 8}] Out[3]= 12.1561 ``` --- Using the product of one-dimensional rules: ```wl In[1]:= NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> {"CartesianRule", Method -> {"GaussKronrodRule", "ClenshawCurtisRule"}}] Out[1]= 328.904 ``` A list of rules is automatically interpreted as a product of rules: ```wl In[2]:= NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> {"GaussKronrodRule", "ClenshawCurtisRule"}] Out[2]= 328.904 ``` Use uniform sampling in ``x`` and nonuniform sampling in ``y``: ```wl In[3]:= ListPlot[Last[Reap[NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> {"NewtonCotesRule", "ClenshawCurtisRule"}, EvaluationMonitor :> Sow[{x, y}], MaxRecursion -> 0 ]]], Frame -> True]//Quiet Out[3]= [image] ``` --- Multivariate integration using a multidimensional symmetric rule: ```wl In[1]:= NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> "MultidimensionalRule"] Out[1]= 328.904 ``` The multidimensional rule uses a sparse symmetric grid of sample points: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> "MultidimensionalRule", EvaluationMonitor :> Sow[{x, y}], MaxRecursion -> 0 ]]], Frame -> True]//Quiet Out[2]= [image] ``` Increase the number of sample points: ```wl In[3]:= ListPlot[Last[Reap[NIntegrate[Log[1 + x ^ 2 + y], {x, 0, 10}, {y, 0, 10}, Method -> {"MultidimensionalRule", "Generators" -> 9}, EvaluationMonitor :> Sow[{x, y}], MaxRecursion -> 0 ]]], Frame -> True]//Quiet Out[3]= [image] ``` Timing can sometimes be improved with a different number of generators: ```wl In[4]:= NIntegrate[Sin[x^2 + y^2], {x, 0, π}, {y, 0, π}]//Timing Out[4]= {0.266, 0.874168} In[5]:= NIntegrate[Sin[x^2 + y^2], {x, 0, π}, {y, 0, π}, Method -> {"MultidimensionalRule", "Generators" -> 9}]//Timing Out[5]= {0.034564, 0.874168} ``` --- Integration of an oscillatory function using a Levin-type collocation rule: ```wl In[1]:= NIntegrate[Sin[(x ^ 3/1 + x)], {x, 0, 100}, Method -> "LevinRule"] Out[1]= 0.696344 ``` Multivariate Levin-type rule: ```wl In[2]:= NIntegrate[Sin[(1/x)]Cos[1000y], {x, 0, 1}, {y, 0, 1}, Method -> "LevinRule"] Out[2]= 0.000416803 ``` ##### Integration Strategies (6) Globally adaptive integration strategy: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "GlobalAdaptive"] Out[1]= 12.1561 ``` Subdivide regions based on largest error until global error is sufficiently small: ```wl In[1]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "GlobalAdaptive", EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0] Out[1]= [image] ``` --- Locally adaptive integration strategy: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "LocalAdaptive"] Out[1]= 12.1561 ``` Subdivide every region until all local errors are sufficiently small: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "LocalAdaptive", EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, AxesOrigin -> {0, 0}] Out[2]= [image] ``` --- Trapezoidal strategy that samples uniformly at increasing density: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "Trapezoidal"] Out[1]= 12.1561 ``` Subdivide entire region and use a lower-order method: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "Trapezoidal", EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, AxesOrigin -> {0, 0}] Out[2]= [image] ``` --- Double-exponential ("tanh-sinh") strategy that samples densely near the endpoints: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "DoubleExponential"] Out[1]= 12.1561 ``` Subdivide entire region, after transformation: ```wl In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "DoubleExponential", EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, AxesOrigin -> {0, 0}] Out[2]= [image] ``` --- Monte Carlo integration with uniformly random sampling points: ```wl In[1]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "MonteCarlo"] Out[1]= 12.1284 In[2]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "MonteCarlo", MaxPoints -> 100, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, AxesOrigin -> {0, 0}]//Quiet Out[2]= [image] ``` With deterministic sequences of sampling points: ```wl In[3]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "QuasiMonteCarlo"] Out[3]= 12.1604 In[4]:= ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "QuasiMonteCarlo", MaxPoints -> 100, EvaluationMonitor :> Sow[{x, Exp[Cos[x]]}]]]], Filling -> 0, AxesOrigin -> {0, 0}]//Quiet Out[4]= [image] ``` Globally adaptive versions of Monte Carlo and quasi Monte Carlo: ```wl In[5]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "AdaptiveMonteCarlo"] Out[5]= 12.0194 In[6]:= NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> "AdaptiveQuasiMonteCarlo"] Out[6]= 12.156 ``` --- Plot sampling points used by different strategies: ```wl In[1]:= Table[ListPlot[Last[Reap[NIntegrate[Exp[Cos[x]], {x, 0, 10}, Method -> str, EvaluationMonitor :> Sow[x]]]], PlotLabel -> str], {str, {"GlobalAdaptive", "LocalAdaptive", "Trapezoidal", "DoubleExponential"}}] Out[1]= {[image], [image], [image], [image]} ``` ##### Symbolic Processing (5) By default some symbolic processing may be performed: ```wl In[1]:= NIntegrate[Exp[-x], {x, 0, 5}] Out[1]= 0.993262 ``` Use the automatic numeric methods, but no symbolic processing: ```wl In[4]:= NIntegrate[Exp[-x], {x, 0, 5}, Method -> {Automatic, "SymbolicProcessing" -> False}] Out[4]= 0.993262 ``` Use an explicit time limit for symbolic processing: ```wl In[5]:= NIntegrate[Exp[-x], {x, 0, 5}, Method -> {Automatic, "SymbolicProcessing" -> 1}] Out[5]= 0.993262 ``` --- Control automatic subdivision of piecewise functions: ```wl In[1]:= NIntegrate[x Ceiling[x], {x, 0, 3}, Method -> {"SymbolicPreprocessing", "SymbolicPiecewiseSubdivision" -> True}] Out[1]= 11. In[2]:= NIntegrate[x Ceiling[x], {x, 0, 3}, Method -> {"SymbolicPreprocessing", "SymbolicPiecewiseSubdivision" -> False}] ``` NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {2.00386}. NIntegrate obtained 11.000949319014198\` and 0.0015217342374763663\` for the integral and error estimates. ```wl Out[2]= 11.0009 ``` Automatic subdivision of piecewise functions usually leads to fewer function evaluations: ```wl In[3]:= Table[ListPlot[Last[Reap[NIntegrate[x Ceiling[x], {x, 0, 3}, Method -> {"SymbolicPreprocessing", o}, EvaluationMonitor :> Sow[{x, x Ceiling[x]}]]]], PlotLabel -> Style[o, Small], Filling -> Axis], {o, {"SymbolicPiecewiseSubdivision" -> True, "SymbolicPiecewiseSubdivision" -> False}}]//Quiet Out[3]= {[image], [image]} ``` --- Control automatic simplification of even and odd integrands: ```wl In[1]:= NIntegrate[x ^ 2 + Cos[x], {x, -3, 3}, Method -> {"SymbolicPreprocessing", "EvenOddSubdivision" -> True}] Out[1]= 18.2822 In[2]:= NIntegrate[x ^ 2 + Cos[x], {x, -3, 3}, Method -> {"SymbolicPreprocessing", "EvenOddSubdivision" -> False}] Out[2]= 18.2822 ``` Automatic simplification usually leads to fewer function evaluations: ```wl In[3]:= Table[ListPlot[Last[Reap[NIntegrate[x ^ 2 + Cos[x], {x, -3, 3}, Method -> {"SymbolicPreprocessing", o}, EvaluationMonitor :> Sow[{x, x ^ 2 + Cos[2x]}]]]], PlotLabel -> o, Filling -> Axis], {o, {"EvenOddSubdivision" -> True, "EvenOddSubdivision" -> False}}] Out[3]= {[image], [image]} ``` --- Control automatic method selection for highly oscillatory functions: ```wl In[1]:= NIntegrate[Sin[x], {x, 0, 10}, Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> True}] Out[1]= 1.83907 In[2]:= NIntegrate[Sin[x], {x, 0, 10}, Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}] Out[2]= 1.83907 ``` Specialized methods lead to fewer function evaluations for highly oscillatory functions: ```wl In[3]:= Table[ListPlot[Last[Reap[NIntegrate[Sin[x], {x, 0, 50}, Method -> {"SymbolicPreprocessing", o}, EvaluationMonitor :> Sow[{x, Sin[x]}]]]], PlotLabel -> o, Filling -> Axis], {o, {"OscillatorySelection" -> True, "OscillatorySelection" -> False}}] Out[3]= {[image], [image]} ``` Switching detection off can save time for non-oscillatory functions: ```wl In[4]:= Timing[Do[NIntegrate[Exp[x], {x, 0, 1}, Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> True}], {10}]] Out[4]= {0.296, Null} In[5]:= Timing[Do[NIntegrate[Exp[x], {x, 0, 1}, Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}], {10}]] Out[5]= {0.172, Null} ``` --- Control automatic subdivision at nodes of interpolating functions: ```wl In[1]:= f = Interpolation[{1, 2, 4, 5, 2, 3}]; In[2]:= NIntegrate[f[x] ^ 2, {x, 0, 5}, Method -> {"SymbolicPreprocessing", "InterpolationPointsSubdivision" -> True}] Out[2]= 49.6771 In[3]:= NIntegrate[f[x] ^ 2, {x, 0, 5}, Method -> {"SymbolicPreprocessing", "InterpolationPointsSubdivision" -> False}]//Quiet Out[3]= 49.6771 ``` Subdivision of interpolating functions can lead to fewer evaluations: ```wl In[4]:= Table[ListPlot[Last[Reap[NIntegrate[f[x] ^ 2, {x, 0, 5}, Method -> {"SymbolicPreprocessing", o}, EvaluationMonitor :> Sow[{x, f[x] ^ 2}]]]], PlotLabel -> Style[o, Small], Filling -> Axis], {o, {"InterpolationPointsSubdivision" -> True, "InterpolationPointsSubdivision" -> False}}]//Quiet Out[4]= {[image], [image]} ``` Subdivision may be unnecessary for an accurate interpolation of a smooth function: ```wl In[5]:= g = FunctionInterpolation[Sin[x], {x, 0, 5}]; In[6]:= Table[ListPlot[Last[Reap[NIntegrate[g[x] ^ 2, {x, 0, 5}, Method -> {"SymbolicPreprocessing", o}, EvaluationMonitor :> Sow[{x, g[x] ^ 2}]]]], PlotLabel -> Style[o, Small], Filling -> Axis], {o, {"InterpolationPointsSubdivision" -> True, "InterpolationPointsSubdivision" -> False}}]//Quiet Out[6]= {[image], [image]} ``` #### MinRecursion (1) ``NIntegrate`` may miss sharp peaks of integrands: ```wl In[1]:= NIntegrate[Exp[-100(x^2 + y^2)], {x, -50, 60}, {y, -50, 60}] ``` NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in x near {x,y} = {23.3334,42.8536}. NIntegrate obtained 0.\` and 0.\` for the integral and error estimates. ```wl Out[1]= 0. ``` Increasing ``MinRecursion`` forces a finer subdivision of the integration region: ```wl In[2]:= NIntegrate[Exp[-100(x^2 + y^2)], {x, -50, 60}, {y, -50, 60}, MinRecursion -> 4] Out[2]= 0.0314159 ``` #### PrecisionGoal (1) The number of samples used to evaluate $\int _0^1\int _0^1\frac{1}{\sqrt{x+y}}dydx$ for different relative tolerances: ```wl In[1]:= sampleCount[pg_] := Module[{k = 0}, NIntegrate[1 / Sqrt[x + y], {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> pg, EvaluationMonitor :> k++]; k ]; ``` The number of samples needed typically increases exponentially with the ``PrecisionGoal`` : ```wl In[2]:= Table[sampleCount[pg], {pg, 2, 12}] Out[2]= {119, 289, 456, 898, 2122, 5281, 12710, 32430, 82342, 201512, 440957} In[3]:= ListLogPlot[%, Joined -> True] Out[3]= [image] ``` #### WorkingPrecision (1) ``NIntegrate`` can compute integrals using higher working precision: ```wl In[1]:= NIntegrate[Exp[-t ^ 2], {t, 0, 12}, WorkingPrecision -> 100]2 / Sqrt[Pi] Out[1]= 0.999999999999999999999999999999999999999999999999999999999999999864373883079409578721969384340958243 ``` The ``PrecisionGoal`` used is 10 less than the ``WorkingPrecision``: ```wl In[2]:= NIntegrate[Exp[-t ^ 2], {t, 0, 12}, WorkingPrecision -> 100, PrecisionGoal -> 90]2 / Sqrt[Pi] Out[2]= 0.999999999999999999999999999999999999999999999999999999999999999864373883079409578721969384340958243 ``` ### Applications (26) #### Basic Applications (2) Compute an integral that has no closed-form solution: ```wl In[1]:= Integrate[(BesselJ[0, x]/1 + x), {x, 0, 1}] Out[1]= Subsuperscript[∫, 0, 1](BesselJ[0, x]/1 + x)\[DifferentialD]x In[2]:= NIntegrate[(BesselJ[0, x]/1 + x), {x, 0, 1}] Out[2]= 0.646543 ``` --- Integrate a discrete set of data using ``Interpolation`` : ```wl In[1]:= data = Table[{x, Sin[(1/x + 1 / 2)]}, {x, 0, 10, 0.2}]; In[2]:= NIntegrate[Interpolation[data][x], {x, 0, 10}] Out[2]= 2.74223 ``` ``Integrate`` also works on interpolating functions: ```wl In[3]:= Integrate[Interpolation[data][x], {x, 0, 10}] Out[3]= 2.74223 ``` Plot the data and the interpolation: ```wl In[4]:= {ListPlot[data, Filling -> 0, AxesOrigin -> {0, 0}], Plot[Interpolation[data][x], {x, 0, 10}, Filling -> 0, AxesOrigin -> {0, 0}]} Out[4]= {[image], [image]} ``` #### Probability and Expectation (3) Compute the probability that $(x-2)^2<1$ when $x$ follows a standard normal distribution: ```wl In[1]:= pdf = PDF[NormalDistribution[0, 1], x] Out[1]= (E^-(x^2/2)/Sqrt[2 π]) In[2]:= NIntegrate[Boole[(x - 2) ^ 2 < 1]pdf, {x, -∞, ∞}] Out[2]= 0.157305 ``` Use ``NProbability`` directly: ```wl In[3]:= NProbability[(x - 2) ^ 2 < 1, x\[Distributed]NormalDistribution[0, 1]] Out[3]= 0.157305 ``` --- Compute the expectation of $\sqrt{| x| }$ when $x$ follows a standard Cauchy distribution: ```wl In[1]:= pdf = PDF[CauchyDistribution[0, 1], x] Out[1]= (1/π (1 + x^2)) In[2]:= NIntegrate[Sqrt[Abs[x]] pdf, {x, -∞, ∞}] Out[2]= 1.41421 ``` Use ``NExpectation`` directly: ```wl In[3]:= NExpectation[Sqrt[Abs[x]], x\[Distributed]CauchyDistribution[0, 1]] Out[3]= 1.41421 ``` --- Compute the cumulative distribution function (CDF) from the probability density function (PDF): ```wl In[1]:= pdf = PDF[NormalDistribution[0, 1], ξ] Out[1]= (E^-(ξ^2/2)/Sqrt[2 π]) In[2]:= cdf[x_ ? NumberQ] := NIntegrate[pdf, {ξ, -∞, x}] In[3]:= Plot[cdf[x], {x, -3, 3}] Out[3]= [image] ``` Compute a quantile value by solving an equation: ```wl In[4]:= quantile[p_ ? NumberQ] := x /. FindRoot[cdf[x] == p, {x, 0}] In[5]:= quantile[0.2] Out[5]= -0.841621 In[6]:= Plot[quantile[p], {p, 0.05, 0.95}, PlotRange -> {{0, 1}, Automatic}] Out[6]= [image] ``` Use ``Quantile`` directly: ```wl In[7]:= Quantile[NormalDistribution[0, 1], 0.2] Out[7]= -0.841621 ``` #### Lengths, Areas, and Volumes (7) Compute the area between two curves as a one-dimensional integral of their difference: ```wl In[1]:= NIntegrate[Exp[Sin[x]] - Sin[Sin[x]], {x, 0, 2}] Out[1]= 2.98948 ``` As a two-dimensional integral over the region bounded by them: ```wl In[2]:= NIntegrate[1, {x, 0, 2}, {y, Sin[Sin[x]], Exp[Sin[x]]}] Out[2]= 2.98948 ``` Plot the computed area: ```wl In[3]:= Plot[{Exp[Sin[x]], Sin[Sin[x]]}, {x, 0, 2}, Filling -> 1 -> {2}] Out[3]= [image] ``` --- Compute the area of a disk that is given implicitly: ```wl In[1]:= NIntegrate[Boole[x ^ 2 + y ^ 2 ≤ 1], {x, -1, 1}, {y, -1, 1}] Out[1]= 3.14159 In[2]:= RegionPlot[x ^ 2 + y ^ 2 ≤ 1, {x, -1, 1}, {y, -1, 1}] Out[2]= [image] ``` Disk annulus: ```wl In[3]:= NIntegrate[Boole[1 / 4 ≤ x ^ 2 + y ^ 2 ≤ 1], {x, -1, 1}, {y, -1, 1}] Out[3]= 2.35619 In[4]:= RegionPlot[1 / 4 ≤ x ^ 2 + y ^ 2 ≤ 1, {x, -1, 1}, {y, -1, 1}] Out[4]= [image] ``` Ellipse: ```wl In[5]:= NIntegrate[Boole[x ^ 2 + (2y) ^ 2 ≤ 1], {x, -1, 1}, {y, -1, 1}] Out[5]= 1.5708 In[6]:= RegionPlot[x ^ 2 + (2y) ^ 2 ≤ 1, {x, -1, 1}, {y, -1, 1}] Out[6]= [image] ``` Ellipse annulus: ```wl In[7]:= NIntegrate[Boole[1 / 4 ≤ x ^ 2 + (2y) ^ 2 ≤ 1], {x, -1, 1}, {y, -1, 1}] Out[7]= 1.1781 In[8]:= RegionPlot[1 / 4 ≤ x ^ 2 + (2y) ^ 2 ≤ 1, {x, -1, 1}, {y, -1, 1}] Out[8]= [image] ``` Disk segment: ```wl In[9]:= NIntegrate[Boole[x ^ 2 + y ^ 2 ≤ 1 && y ≥ 0 && y ≥ -x], {x, -1, 1}, {y, -1, 1}] Out[9]= 1.1781 In[10]:= RegionPlot[x ^ 2 + y ^ 2 ≤ 1 && y ≥ 0 && y ≥ -x, {x, -1, 1}, {y, -1, 1}] Out[10]= [image] ``` --- Simple regions including a ball: ```wl In[1]:= NIntegrate[Boole[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[1]= 21.7656 In[2]:= RegionPlot3D[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None] Out[2]= [image] ``` Half of a spherical shell: ```wl In[3]:= NIntegrate[Boole[1 ≤ x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[3]= 8.7884 + 4.585245263883044`*^-59 I In[4]:= RegionPlot3D[1 ≤ x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None, PlotPoints -> 50] Out[4]= [image] ``` Ellipsoid: ```wl In[5]:= NIntegrate[Boole[x ^ 2 + y ^ 2 + 1 / 2z ^ 2 ≤ 3], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}] Out[5]= 30.7812 In[6]:= RegionPlot3D[x ^ 2 + y ^ 2 + 1 / 2z ^ 2 ≤ 3, {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, Mesh -> None, PlotPoints -> 35] Out[6]= [image] ``` Half of an ellipsoidal shell: ```wl In[7]:= NIntegrate[Boole[1 ≤ x ^ 2 + y ^ 2 + 1 / 2z ^ 2 ≤ 3 && y ≥ 0], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}] Out[7]= 12.4287 + 6.484516038990402`*^-59 I In[8]:= RegionPlot3D[1 ≤ x ^ 2 + y ^ 2 + 1 / 2z ^ 2 ≤ 3 && y ≥ 0, {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, Mesh -> None, PlotPoints -> 50] Out[8]= [image] ``` Spherical wedge: ```wl In[9]:= NIntegrate[Boole[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0 && y ≥ x + z], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}] Out[9]= 7.57348 + 2.0875171654287265`*^-59 I In[10]:= RegionPlot3D[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0 && y ≥ x + z, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None, PlotPoints -> 50] Out[10]= [image] ``` --- Compute the area of a region defined by a closed parametric curve using Green's theorem: ```wl In[1]:= {x, y} = {Cos[u], Sin[u]}; In[2]:= NIntegrate[x D[y, u], {u, 0, 2π}] Out[2]= 3.14159 In[3]:= ParametricPlot[{x, y}, {u, 0, 2π}] Out[3]= [image] ``` Ellipse: ```wl In[4]:= {x, y} = {2Cos[u], Sin[u]}; In[5]:= NIntegrate[x D[y, u], {u, 0, 2π}] Out[5]= 6.28319 In[6]:= ParametricPlot[{x, y}, {u, 0, 2Pi}] Out[6]= [image] ``` A periodic boundary curve: ```wl In[7]:= {x, y} = {Cos[u] + 1 / 7Cos[7u + Pi / 3], Sin[u] + 1 / 7Sin[7u]}; In[8]:= NIntegrate[x D[y, u], {u, 0, 2π}] Out[8]= 3.36599 In[9]:= ParametricPlot[{x, y}, {u, 0, 2Pi}] Out[9]= [image] ``` --- Compute the length of a parametric curve: ```wl In[1]:= {x, y} = {Cos[u], Sin[u]}; In[2]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2], {u, 0, 2π}] Out[2]= 6.28319 In[3]:= ParametricPlot[{x, y}, {u, 0, 2π}] Out[3]= [image] ``` Ellipse: ```wl In[4]:= {x, y} = {2Cos[u], Sin[u]}; In[5]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2], {u, 0, 2π}] Out[5]= 9.68845 In[6]:= ParametricPlot[{x, y}, {u, 0, 2π}] Out[6]= [image] ``` A periodic boundary curve: ```wl In[7]:= {x, y} = {Cos[u] + 1 / 7Cos[7u + Pi / 3], Sin[u] + 1 / 7Sin[7u]}; In[8]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2], {u, 0, 2π}] Out[8]= 8.09043 In[9]:= ParametricPlot[{x, y}, {u, 0, 2Pi}] Out[9]= [image] ``` A circle in 3D: ```wl In[10]:= {x, y, z} = {1, 0, 0}Cos[u] + Normalize[{0, 1, 1}]Sin[u]; In[11]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2 + D[z, u] ^ 2], {u, 0, 2π}] Out[11]= 6.28319 In[12]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}] Out[12]= [image] ``` Ellipse in 3D: ```wl In[13]:= {x, y, z} = {1, 0, 0}2Cos[u] + Normalize[{0, 1, 1}]Sin[u]; In[14]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2 + D[z, u] ^ 2], {u, 0, 2π}] Out[14]= 9.68845 In[15]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}] Out[15]= [image] ``` More general curve: ```wl In[16]:= {x, y, z} = {1, 0, 0}(Cos[u] + 1 / 7Cos[7u + Pi / 3]) + Normalize[{1, 1, 1}](Sin[u] + 1 / 7Sin[7u]); In[17]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2 + D[z, u] ^ 2], {u, 0, 2π}] Out[17]= 7.18498 In[18]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}] Out[18]= [image] ``` A toroidal spring curve: ```wl In[19]:= {x, y, z} = {Cos[u] (2 + Cos[8 u]), (2 + Cos[8 u]) Sin[u], Sin[8 u]}; In[20]:= NIntegrate[Sqrt[D[x, u] ^ 2 + D[y, u] ^ 2 + D[z, u] ^ 2], {u, 0, 2π}] Out[20]= 51.9914 In[21]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}] Out[21]= [image] ``` --- Compute the area of a parametrically defined surface: ```wl In[1]:= {x, y, z} = {Cos[u]Sin[v], Sin[u]Sin[v], Cos[v]}; In[2]:= NIntegrate[Norm[Cross[D[{x, y, z}, u], D[{x, y, z}, v]]], {u, 0, 2π}, {v, 0, π}] Out[2]= 12.5664 In[3]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, π}] Out[3]= [image] ``` Ellipsoidal surface: ```wl In[4]:= {x, y, z} = {2Cos[u]Sin[v], Sin[u]Sin[v], Cos[v]}; In[5]:= NIntegrate[Norm[Cross[D[{x, y, z}, u], D[{x, y, z}, v]]], {u, 0, 2π}, {v, 0, π}] Out[5]= 21.4784 In[6]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, π}] Out[6]= [image] ``` Toroidal surface: ```wl In[7]:= {x, y, z} = {(2 + Cos[v])Cos[u], (2 + Cos[v])Sin[u], Sin[v]}; In[8]:= NIntegrate[Norm[Cross[D[{x, y, z}, u], D[{x, y, z}, v]]], {u, 0, 2π}, {v, 0, 2π}] Out[8]= 78.9568 In[9]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}] Out[9]= [image] ``` General parametric surface: ```wl In[10]:= {x, y, z} = {(u/2) + v, u - (v/2), Sin[u v]}; In[11]:= NIntegrate[Norm[Cross[D[{x, y, z}, u], D[{x, y, z}, v]]], {u, 0, 3}, {v, 0, 3}] Out[11]= 19.477 In[12]:= ParametricPlot3D[{x, y, z}, {v, 0, 3}, {u, 0, 3}] Out[12]= [image] ``` --- Compute the volume enclosed by a parametric surface using the divergence theorem: ```wl In[1]:= {x, y, z} = {Cos[u]Sin[v], Sin[u]Sin[v], Cos[v]}; In[2]:= NIntegrate[({x, y, z}/3).Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, π}]//Abs Out[2]= 4.18879 In[3]:= 4 / 3 Pi//N Out[3]= 4.18879 In[4]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}] Out[4]= [image] ``` Ellipsoidal volume: ```wl In[5]:= {x, y, z} = {4Cos[u]Sin[v], 3Sin[u]Sin[v], 2Cos[v]}; In[6]:= NIntegrate[({x, y, z}/3).Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, π}]//Abs Out[6]= 100.531 In[7]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}] Out[7]= [image] ``` Toroidal volume: ```wl In[8]:= {x, y, z} = {(2 + Cos[v])Cos[u], (2 + Cos[v])Sin[u], Sin[v]}; In[9]:= NIntegrate[({x, y, z}/3).Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, 2π}] Out[9]= 39.4784 In[10]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}] Out[10]= [image] ``` Volume of a general parametric surface: ```wl In[11]:= {x, y, z} = {(2 + Cos[v])Cos[u], (2 + Cos[v])Sin[u], Sin[v](5 / 4 - Cos[u])}; In[12]:= NIntegrate[({x, y, z}/3).Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, 2π}] Out[12]= 49.348 In[13]:= ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}] Out[13]= [image] ``` #### Line, Surface, and Volume Integrals (3) Work done by a cyclic thermodynamic process: ```wl In[1]:= {p, v} = {5 + (3 + Sin[3t])Cos[t], 6 + (3 + Sin[t / 2] ^ 2)Sin[t]}; In[2]:= NIntegrate[p D[v, t], {t, 0, 2π}] Out[2]= 32.9867 ``` Visualize cycle on a pressure-volume diagram: ```wl In[3]:= ParametricPlot[{p, v}, {t, 0, 2π}, PlotRange -> {{0, 9}, {0, 10}}, AxesLabel -> {"V", "P"}] /. Line[p___] :> {Arrowheads[Medium], Arrow[p]} Out[3]= [image] ``` --- Mass and center of mass of a region with uniform density: ```wl In[1]:= m = NIntegrate[Boole[2 < Abs[2x] + y ^ 2 < 3 && 0 < y < 3], {x, -∞, ∞}, {y, -∞, ∞}] Out[1]= 1.57848 In[2]:= cm = (1/m)NIntegrate[{x, y}Boole[2 < Abs[2x] + y ^ 2 < 3 && 0 < y < 3], {x, -∞, ∞}, {y, -∞, ∞}] Out[2]= {0., 0.791899} ``` Visualize the center of mass: ```wl In[3]:= RegionPlot[2 < Abs[2x] + y ^ 2 < 3 && 0 < y < 3, {x, -2, 2}, {y, -1, 3}, AxesOrigin -> cm, Axes -> True, AxesStyle -> Dashed] Out[3]= [image] ``` Three-dimensional region: ```wl In[4]:= m = NIntegrate[Boole[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[4]= 10.8828 In[5]:= cm = (1/m)NIntegrate[{x, y, z}Boole[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Out[5]= {0., 0.649519, 0.} ``` Visualize the center of mass: ```wl In[6]:= Show[RegionPlot3D[x ^ 2 + y ^ 2 + z ^ 2 ≤ 3 && y ≥ 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None, PlotPoints -> 50, PlotStyle -> Opacity[1 / 2]], Graphics3D@Line[Transpose[{cm + IdentityMatrix[3] / 3, cm - IdentityMatrix[3] / 3}, {2, 3, 1}]]] Out[6]= [image] ``` --- Center of mass of a region enclosed by a parametric surface using Stokes's theorem: ```wl In[1]:= {x, y, z} = {(2 + Cos[v])Cos[u], (2 + Cos[v])Sin[u], Sin[v](5 / 4 - Cos[u])}; In[2]:= m = NIntegrate[({x, y, z}/3).Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, 2π}] Out[2]= 49.348 In[3]:= cm = (1/m)NIntegrate[DiagonalMatrix[({x, y, z} ^ 2/2)].Cross[D[{x, y, z}, u], D[{x, y, z}, v]], {u, 0, 2π}, {v, 0, 2π}, AccuracyGoal -> 4] Out[3]= {-0.85, -5.316794286643837`*^-16, -3.624976031461681`*^-16} ``` Visualize center of mass: ```wl In[4]:= Show[ParametricPlot3D[{x, y, z}, {u, 0, 2π}, {v, 0, 2π}, PlotStyle -> Opacity[1 / 2], Mesh -> None], Graphics3D@Line[Transpose[{cm + 10IdentityMatrix[3] / 3, cm - 10IdentityMatrix[3]}, {2, 3, 1}]]] Out[4]= [image] ``` #### Integral Transforms (4) Compute integral transforms including Fourier transform: ```wl In[1]:= ft[f_, x_, ω_ ? NumberQ] := (1/Sqrt[2π])NIntegrate[f Exp[I ω x], {x, -∞, ∞}, AccuracyGoal -> 2] In[2]:= ft[Exp[-x ^ 2], x, 0.5] Out[2]= 0.664265 + 0. I In[3]:= Plot[Abs[ft[Exp[-x ^ 2], x, ω]], {ω, -3, 3}] Out[3]= [image] ``` Laplace transform: ```wl In[4]:= lt[f_, x_, s_ ? NumberQ] := NIntegrate[f Exp[-s x], {x, 0, ∞}] In[5]:= lt[x ^ 2, x, 0.5] Out[5]= 16. In[6]:= Plot[Abs@lt[x ^ 2, x, 0.5 + I ω], {ω, -3, 3}] Out[6]= [image] ``` Mellin transform: ```wl In[7]:= mt[f_, x_, s_] := NIntegrate[f x ^ (s - 1), {x, 0, ∞}] In[8]:= mt[Sin[x], x, 0.5] Out[8]= 1.25331 In[9]:= Plot[mt[Sin[x], x, s], {s, -1, 1}] Out[9]= [image] ``` Hilbert transform: ```wl In[10]:= ht[f_, x_, τ_] := NIntegrate[f (1/π (x - τ)), {x, -∞, τ, ∞}, Method -> {"PrincipalValue"}] In[11]:= ht[Sinc[x], x, 0.5] Out[11]= -0.244835 In[12]:= Plot[ht[Sinc[x], x, τ], {τ, -5, 5}] Out[12]= [image] ``` Hartley transform: ```wl In[13]:= ht[f_, x_, ω_] := (1/Sqrt[2π])NIntegrate[f(Cos[ω x] + Sin[ω x]), {x, -∞, ∞}] In[14]:= ht[Exp[-x ^ 2], x, 0.5] Out[14]= 0.664265 In[15]:= Plot[ht[Exp[-x ^ 2], x, ω], {ω, -3, 3}] Out[15]= [image] ``` --- Find the Fourier coefficients of a periodic function on ``[0, 1]`` : ```wl In[1]:= a[i_] := 2NIntegrate[x^2 Cos[2π i x], {x, 0, 1}] In[2]:= b[i_] := 2 NIntegrate[x^2 Sin[2π i x], {x, 0, 1}] In[3]:= {DiscretePlot[a[i], {i, 0, 10}], DiscretePlot[b[i], {i, 1, 10}]} Out[3]= {[image], [image]} ``` Show the approximate inverse transform: ```wl In[4]:= Table[Plot[Evaluate[{x^2, a[0] / 2 + Subsuperscript[∑, i = 1, n]a[i]Cos[2π i x] + Subsuperscript[∑, i = 1, n]b[i]Sin[2π i x]}], {x, 0, 1}, Ticks -> None, PlotLabel -> n], {n, 2, 8, 2}] Out[4]= [image] ``` --- Compute a quadratic fractional Fourier transform: ```wl In[1]:= fract[α_, w_] := Sqrt[(1 - I Cot[α]) / (2Pi)]E ^ (I Cot[α] w ^ 2 / 2) NIntegrate[E ^ (-I Csc[α] w t + I Cot[α] t ^ 2 / 2) UnitBox[t], {t, -Infinity, Infinity}] In[2]:= Table[DiscretePlot[{Re[fract[k Pi, w]], Im[fract[k Pi, w]]}, {w, -2, 2, 0.05}, Joined -> True, PlotStyle -> {Red, Blue}, PlotRange -> All, Filling -> Axis], {k, {0.02, 0.09, 3.2, 20.9}}] Out[2]= [image] ``` --- Fraunhofer integral for amplitude of waves diffracted by a square aperture: ```wl In[1]:= a[x_, y_, z_] := (1 /2π z)Abs[NIntegrate[Exp[(x s + y t/I z)], {s, 0, 1}, {t, 0, 1}, PrecisionGoal -> 2]] ``` Diffraction pattern near optical axis: ```wl In[2]:= MatrixPlot[ParallelTable[a[x, y, 1], {x, -10, 10}, {y, -10, 20}]] Out[2]= [image] ``` #### Integral Representations of Special Functions (3) Integral representation for Bessel function of the first kind $J_n(z)$ on the real line: ```wl In[1]:= j[n_, x_ ? NumberQ] := NIntegrate[(Cos[t n - x Sin[t]]/π), {t, 0, π}] ``` Compare with built-in function ``BesselJ[n, x]`` : ```wl In[2]:= Table[Plot[{f[0, x], f[1, x]}, {x, 0, 10}, PlotLabel -> Row[{f[0, x], f[1, x]}, ","]], {f, {j, BesselJ}}] Out[2]= {[image], [image]} ``` --- Gamma function $\Gamma (z)$ on the right half of the complex plane: ```wl In[1]:= g[z_ ? NumberQ] := NIntegrate[t ^ (z - 1)Exp[-t], {t, 0, ∞}] ``` Compare with built-in function ``Gamma[z]`` : ```wl In[2]:= Table[ContourPlot[Abs[f[x + I y]], {x, 1, 3}, {y, -1, 1}, PlotPoints -> 5, PlotLabel -> f[z]], {f, {g, Gamma}}] Out[2]= {[image], [image]} ``` --- Incomplete elliptic integral of the second kind $E(z|m)$ on the real line: ```wl In[1]:= e[z_ ? NumberQ, m_ ? NumberQ] := NIntegrate[Sqrt[1 - m Sin[t] ^ 2], {t, 0, z}] ``` Compare with built-in function ``EllipticE[z, m]`` : ```wl In[2]:= Table[Plot[{f[z, 1], f[2, z]}, {z, 0, 1}, PlotLabel -> Row[{f[z, 1], f[2, z]}, ","]], {f, {e, EllipticE}}] Out[2]= {[image], [image]} ``` #### Function Norms and Inner Products (2) $L^p$ norm of a function: ```wl In[1]:= lp[p_ ? NumberQ, f_, x_] := NIntegrate[Abs[f] ^ p, {x, -∞, ∞}] ^ (1 / p) In[2]:= Table[lp[p, Exp[-5x ^ 2], x], {p, 5}] Out[2]= {0.792665, 0.748665, 0.770625, 0.793442, 0.812683} ``` Plot $L^p$ norm as a function of $p$ : ```wl In[3]:= Table[Plot[lp[p, f, x], {p, 1, 5}, PlotLabel -> Subscript[Norm[f], p], AxesLabel -> {p}], {f, {Exp[-x ^ 2], Exp[-5x ^ 2], Sin[x]Exp[-5x ^ 2]}}] Out[3]= {[image], [image], [image]} ``` Minimize $\| f\| _p$ as a function of $p$ : ```wl In[4]:= FindMinimum[lp[p, Exp[-5x ^ 2], x], {p, 1, 5}] Out[4]= {0.746209, {p -> 1.70795}} ``` --- $L^p$ inner product with respect to weight function $w$ for functions defined on $a 8] ``` Orthogonality of Legendre polynomials $P_n(x)$ on $-1 c + Exp[I t], {t, 0, 2π}] In[2]:= res[(2/z - 1), {z, 1}] Out[2]= 2. + 0. I ``` Compare with exact value: ```wl In[3]:= res[(Sin[z]/(z - 2) ^ 3), {z, 2}] Out[3]= -0.454649 - 1.2090443197003527`*^-16 I In[4]:= Residue[(Sin[z]/(z - 2) ^ 3), {z, 2}]//N Out[4]= -0.454649 ``` --- Numerical derivative of a function at $z=a$ as an integral over a contour enclosing $a$ : ```wl In[1]:= d[f_, z_, a_] := (1/2π I)NIntegrate[(f/(z - a) ^ 2) Dt[z, t] /. z -> a + Exp[I t], {t, 0, 2π}] In[2]:= {d[Sin[z + Sin[z]], z, 1.5 + I], D[Sin[z + Sin[z]], z] /. z -> 1.5 + I} Out[2]= {-1.97302 + 1.77146 I, -1.97302 + 1.77146 I} ``` Differentiate a function defined numerically: ```wl In[3]:= f[x_ ? NumericQ] := y[1] /. First@NDSolve[{y'[s] == 1 - x y[s] ^ 3, y[0] == x}, y, {s, 1}] In[4]:= Plot[f[x], {x, 0, 4}] Out[4]= [image] In[5]:= d[f[x], x, 2] Out[5]= -0.13696 - 2.208718528794109`*^-17 I ``` ### Properties & Relations (9) ``NIntegrate`` computes a surface integral for lower-dimensional regions: ```wl In[1]:= NIntegrate[1, {x, y}∈Circle[]] Out[1]= 6.28319 In[2]:= NIntegrate[1, {x, y, z}∈Sphere[]] Out[2]= 12.5664 ``` The regions are lower dimensional: ```wl In[3]:= RegionDimension[Circle[]] < RegionEmbeddingDimension[Circle[]] Out[3]= True In[4]:= RegionDimension[Sphere[]] < RegionEmbeddingDimension[Sphere[]] Out[4]= True ``` --- The surface (curve etc.) integral $\int _{x\in S}g(x)$ for a parametric surface $\left\{f_1(\theta ),\ldots ,f_n(\theta )\right\}$ is given by the ordinary integral $\int _{\theta \in \mathcal{R}}G(\theta ) g(f(\theta ))$ where $G(\theta )$ is the Gram determinant: ```wl In[1]:= GramDet[f_, v_] := Sqrt[Det[D[f, {v}]\[Transpose].D[f, {v}]]]; ``` Integrating over a parametric circle curve $\{\cos (\phi ),\sin (\phi )\}$ : ```wl In[2]:= G = Simplify@GramDet[{Cos[ϕ], Sin[ϕ]}, {ϕ}] Out[2]= 1 In[3]:= {NIntegrate[G, {ϕ, 0, 2π}], NIntegrate[1, {x, y}∈Circle[]]} Out[3]= {6.28319, 6.28319} ``` Integrating over a parametric sphere surface $\{\sin (\theta ) \cos (\phi ),\sin (\theta ) \sin (\phi ),\cos (\theta )\}$ : ```wl In[4]:= G = Simplify@GramDet[{Cos[ϕ]Sin[θ], Sin[ϕ]Sin[θ], Cos[θ]}, {ϕ, θ}] Out[4]= Sqrt[Sin[θ]^2] In[5]:= {NIntegrate[G, {ϕ, 0, 2π}, {θ, 0, π}], NIntegrate[1, {x, y, z}∈Sphere[]]} Out[5]= {12.5664, 12.5664} ``` --- Integrating over a zero-dimensional region corresponds to summing over the points: ```wl In[1]:= NIntegrate[x y, {x, y}∈Point[{{1, 2}, {3, 1}}]] Out[1]= 5. In[2]:= Sum[p[[1]]p[[2]], {p, {{1, 2}, {3, 1}}}] Out[2]= 5 ``` --- When a closed-form solution is available, ``Integrate`` can be used instead of ``NIntegrate`` : ```wl In[1]:= exact = Integrate[(1/Sqrt[x + y]), {x, 0, 1}, {y, 0, 1}] Out[1]= (8/3) (-1 + Sqrt[2]) ``` The result of ``NIntegrate`` is close to the exact value: ```wl In[2]:= exact - NIntegrate[(1/Sqrt[x + y]), {x, 0, 1}, {y, 0, 1}] Out[2]= -4.57992288538378`*^-9 ``` --- ``NDSolve`` can be used instead of ``NIntegrate`` : ```wl In[1]:= sol = NDSolve[{D[f[x], x] == Sqrt[x], f[0] == 0}, f, {x, 0, 1}][[1, 1]] Out[1]= f -> InterpolatingFunction[{{0., 1.}}, {5, 7, 1, {96}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«1114»"], {Developer`PackedArrayForm, {0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48 ... 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186, 188, 190, 192}, CompressedData["«2135»"]}, {Automatic}] In[2]:= 2 / 3 - (f /. sol)[1] Out[2]= -3.640009060834615`*^-8 ``` ``NIntegrate`` will typically give a more precise result since it uses global error control: ```wl In[3]:= 2 / 3 - NIntegrate[Sqrt[x], {x, 0, 1}] Out[3]= -3.4416913763379853`*^-15 ``` --- Approximate the tail of a sum using the Euler–Maclaurin formula: ```wl In[1]:= f[x_] := 1 / x ^ 2; In[2]:= Sum[f[x], {x, 1, 15}] + NIntegrate[f[x], {x, 15, ∞}] - (f[15] + f[∞]) / 2 + Sum[(BernoulliB[2i](f^(2i - 1)[∞] - f^(2i - 1)[15])) / (2i)!, {i, 1, 6}] Out[2]= 1.64493 ``` The Euler–Maclaurin formula approximation compares well with the exact result: ```wl In[3]:= Sum[1 / x ^ 2, {x, 1, ∞}] - % Out[3]= -6.661338147750939`*^-16 ``` The Euler–Maclaurin method of ``NSum`` uses ``Integrate`` or ``NIntegrate`` : ```wl In[4]:= NSum[1 / x^2, {x, 1, ∞}, Method -> {"EulerMaclaurin", Method -> NIntegrate}] Out[4]= 1.64493 ``` --- Find a cumulative distribution function from a probability density function: ```wl In[1]:= NIntegrate[PDF[NormalDistribution[2, 5], y], {y, -∞, 4.}] Out[1]= 0.655422 ``` ``CDF`` gives the same result: ```wl In[2]:= CDF[NormalDistribution[2, 5], 4.] Out[2]= 0.655422 ``` --- The ``Residue`` of a function at a pole is equivalent to a contour integral: ```wl In[1]:= Residue[x ^ 2 / (x ^ 2 + 1) ^ 2, {x, I}] Out[1]= -(I/4) ``` Numerically integrate on a contour that includes the real axis from ``-∞`` to ``∞`` : ```wl In[2]:= (1/2 π I) NIntegrate[x ^ 2 / (x ^ 2 + 1) ^ 2, {x, -∞, ∞}] Out[2]= -0.25 I ``` --- Several special functions are defined as integrals: ```wl In[1]:= Block[{x = 10.}, {NIntegrate[1 / t, {t, 1, x}], (2/Sqrt[π])NIntegrate[Exp[-t ^ 2], {t, 0, x}], NIntegrate[t ^ (x - 1) / E ^ t, {t, 0, Infinity}]}] Out[1]= {2.30259, 1., 362880.} In[2]:= Block[{x = 10.}, {Log[x], Erf[x], Gamma[x]}] Out[2]= {2.30259, 1., 362880.} ``` For high precision, using the special function is usually preferable: ```wl In[3]:= N[Erf[12], 100] Out[3]= 0.999999999999999999999999999999999999999999999999999999999999999864373883079409578721969384340958243 In[4]:= (2/Sqrt[π])NIntegrate[Exp[-t ^ 2], {t, 0, 12}, PrecisionGoal -> 100, WorkingPrecision -> 110] Out[4]= 0.9999999999999999999999999999999999999999999999999999999999999998643738830794095787219693843409582427333217767 ``` ### Possible Issues (6) For integrands with complicated variation, many levels of adaptive recursion may be required: ```wl In[1]:= NIntegrate[(Cos[20000x]/Sqrt[x]), {x, 0, 1}] ``` NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {0.}. NIntegrate obtained 0.00579916354919178\` and 0.014875547652988541\` for the integral and error estimates. ```wl Out[1]= 0. ``` Specifying higher recursive subdivision levels gives a convergent result: ```wl In[2]:= NIntegrate[(Cos[20000x]/Sqrt[x]), {x, 0, 1}, MaxRecursion -> 20] Out[2]= 0.00889137 ``` --- One-dimensional integral with a singularity inside the integration region: ```wl In[1]:= NIntegrate[Log[(1 - x) ^ 2], {x, 0, 2}] ``` NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {0.999969}. NIntegrate obtained -3.99995 and 0.0002785157244333794\` for the integral and error estimates. ```wl Out[1]= -3.99995 ``` Specify the location of the singularity to handle it with appropriate transformations: ```wl In[2]:= NIntegrate[Log[(1 - x) ^ 2], {x, 0, 1, 2}] Out[2]= -4. In[3]:= NIntegrate[Log[(1 - x) ^ 2], {x, 0, 2}, Exclusions -> {1}] Out[3]= -4. ``` --- Since the default is to use relative tolerance, integrals that are zero may be computed slowly: ```wl In[1]:= NIntegrate[Sin[x]Cos[y], {x, 0, 2π}, {y, 0, 2π}]//Timing ``` NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained -1.46801\*10^-16 and 1.4974799404526401\`\*^-12 for the integral and error estimates. ```wl Out[1]= {32.397, 0``11.824638986507319} ``` Specifying an absolute tolerance with ``AccuracyGoal`` will reduce the amount of work: ```wl In[2]:= NIntegrate[Sin[x]Cos[y], {x, 0, 2π}, {y, 0, 2π}, AccuracyGoal -> 8]//Timing Out[2]= {0.01, 0``14.058950510203026} ``` --- Higher-dimensional cubature-based integration to default precision may be slow: ```wl In[1]:= NIntegrate[1 / Sqrt[x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9], {x1, 0, 1}, {x2, 0, 1}, {x3, 0, 1}, {x4, 0, 1}, {x5, 0, 1}, {x6, 0, 1}, {x7, 0, 1}, {x8, 0, 1}, {x9, 0, 1}]//Timing Out[1]= {109.507, 0.478548} ``` Decreasing the ``PrecisionGoal`` speeds up the computation: ```wl In[2]:= NIntegrate[1 / Sqrt[x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9], {x1, 0, 1}, {x2, 0, 1}, {x3, 0, 1}, {x4, 0, 1}, {x5, 0, 1}, {x6, 0, 1}, {x7, 0, 1}, {x8, 0, 1}, {x9, 0, 1}, PrecisionGoal -> 2]//Timing Out[2]= {0.01, 0.478554} ``` Using a quasi Monte Carlo method for the same low ``PrecisionGoal`` may be even faster: ```wl In[3]:= NIntegrate[1 / Sqrt[x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9], {x1, 0, 1}, {x2, 0, 1}, {x3, 0, 1}, {x4, 0, 1}, {x5, 0, 1}, {x6, 0, 1}, {x7, 0, 1}, {x8, 0, 1}, {x9, 0, 1}, Method -> "QuasiMonteCarlo"]//Timing Out[3]= {0.01, 0.477587} ``` ``NIntegrate`` automatically uses Monte Carlo integration for dimensions higher than 15: ```wl In[4]:= NIntegrate[1 / (√(x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16)), {x1, 0, 1}, {x2, 0, 1}, {x3, 0, 1}, {x4, 0, 1}, {x5, 0, 1}, {x6, 0, 1}, {x7, 0, 1}, {x8, 0, 1}, {x9, 0, 1}, {x10, 0, 1}, {x11, 0, 1}, {x12, 0, 1}, {x13, 0, 1}, {x14, 0, 1}, {x15, 0, 1}, {x16, 0, 1}]//Timing Out[4]= {0.01, 0.357579} ``` --- It can be time-consuming to compute functions symbolically: ```wl In[1]:= f[x_] := Nest[Sin[# + Sin[2#]]&, x, 20] In[2]:= NIntegrate[f[x], {x, 0, 1}]//Timing Out[2]= {1.20379, 0.947747} ``` Restricting the function definition avoids symbolic evaluation: ```wl In[3]:= g[x_ ? NumericQ] := Nest[Sin[# + Sin[2#]]&, x, 20] In[4]:= NIntegrate[g[x], {x, 0, 1}]//Timing Out[4]= {0.003489, 0.947747} ``` --- Define a function that does numerical integration for a given parameter: ```wl In[1]:= f1[a_] := NIntegrate[(x - a) ^ 2, {x, 0, 10}] ``` Compute with a parameter value of 2: ```wl In[2]:= f1[2] Out[2]= 173.333 ``` Applying the function to a symbolic parameter generates a message from ``NIntegrate`` : ```wl In[3]:= f1[q] ``` NIntegrate::inumr: The integrand (-q+x)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,10}}. ```wl Out[3]= NIntegrate[(x - q)^2, {x, 0, 10}] ``` This can also lead to warnings when the function is used with other numerical functions like ``NMinimize`` : ```wl In[4]:= NMinimize[f1[a], a] ``` NIntegrate::inumr: The integrand (-a+x)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,10}}. NIntegrate::inumr: The integrand (-a+x)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,10}}. NIntegrate::inumr: The integrand (-a+x)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,10}}. General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation. ```wl Out[4]= {83.3333, {a -> 5.}} ``` Define a function that only evaluates when its argument is a numerical value to avoid these messages: ```wl In[5]:= f2[a_ ? NumericQ] := NIntegrate[(x - a) ^ 2, {x, 0, 10}] ``` Compute with a numerical value: ```wl In[6]:= f2[π] Out[6]= 117.87 ``` The function does not evaluate when its argument is non-numerical: ```wl In[7]:= f2[q] Out[7]= f2[q] ``` The function can now be used with other numerical functions such as ``NMinimize`` : ```wl In[8]:= NMinimize[f2[a], a] Out[8]= {83.3333, {a -> 5.}} ``` ### Neat Examples (2) Sampling points for a three-dimensional integration, using a low-order method: ```wl In[1]:= points = Reap[NIntegrate[Boole[(x1)^2 + (x2 - 1)^2 < 1] + Boole[(x1)^2 + (x3 - 1)^2 < 1], {x1, 0, 1}, {x2, 0, 2}, {x3, 0, 2}, Method -> {"GlobalAdaptive", Method -> {"TrapezoidalRule", "Points" -> 3}, "SingularityHandler" -> None}, PrecisionGoal -> 4, EvaluationMonitor :> Sow[{x1, x2, x3}]]][[2, 1]]; In[2]:= Graphics3D[Point[points]] Out[2]= [image] ``` --- Integrate a function with several singular curves: ```wl In[1]:= NIntegrate[Sqrt[(1/Cos[5x]Sin[5y] - 1 / 2)], {x, 0, 5}, {y, 0, 5}, Exclusions -> Cos[5x]Sin[5y] == 1 / 2]//Quiet Out[1]= 6.63973 + 32.5811 I In[2]:= Plot3D[{Re[Sqrt[(1/Cos[5x]Sin[5y] - 1 / 2)]], Im[Sqrt[(1/Cos[5x]Sin[5y] - 1 / 2)]]}, {x, 0, 5}, {y, 0, 5}, Exclusions -> Cos[5x]Sin[5y] == 1 / 2, Mesh -> None, PlotStyle -> {Yellow, Magenta}, PlotPoints -> 35] Out[2]= [image] ``` ## See Also * [`NDSolve`](https://reference.wolfram.com/language/ref/NDSolve.en.md) * [`NSum`](https://reference.wolfram.com/language/ref/NSum.en.md) * [`Integrate`](https://reference.wolfram.com/language/ref/Integrate.en.md) * [`AsymptoticIntegrate`](https://reference.wolfram.com/language/ref/AsymptoticIntegrate.en.md) * [`NLineIntegrate`](https://reference.wolfram.com/language/ref/NLineIntegrate.en.md) * [`NSurfaceIntegrate`](https://reference.wolfram.com/language/ref/NSurfaceIntegrate.en.md) * [`NContourIntegrate`](https://reference.wolfram.com/language/ref/NContourIntegrate.en.md) * [`NExpectation`](https://reference.wolfram.com/language/ref/NExpectation.en.md) * [`NProbability`](https://reference.wolfram.com/language/ref/NProbability.en.md) * [`ArcLength`](https://reference.wolfram.com/language/ref/ArcLength.en.md) * [`Area`](https://reference.wolfram.com/language/ref/Area.en.md) * [`Volume`](https://reference.wolfram.com/language/ref/Volume.en.md) * [`MomentOfInertia`](https://reference.wolfram.com/language/ref/MomentOfInertia.en.md) ## Tech Notes * [Advanced Numerical Integration in the Wolfram Language](https://reference.wolfram.com/language/tutorial/NIntegrateOverview.en.md) * [Numerical Mathematics: Basic Operations](https://reference.wolfram.com/language/tutorial/NumericalCalculations.en.md) * [Numerical Integration](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#12104) * [Numerical Sums, Products, and Integrals](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#3848) * [Numerical Mathematics in the Wolfram Language](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#20411) * [The Uncertainties of Numerical Mathematics](https://reference.wolfram.com/language/tutorial/NumericalOperationsOnFunctions.en.md#14523) * [Implementation notes: Numerical and Related Functions](https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.en.md#23656) ## Related Guides * [`Calculus`](https://reference.wolfram.com/language/guide/Calculus.en.md) * [Geometric Computation](https://reference.wolfram.com/language/guide/GeometricComputation.en.md) * [Functions of Complex Variables](https://reference.wolfram.com/language/guide/FunctionsOfComplexVariables.en.md) * [Solvers over Regions](https://reference.wolfram.com/language/guide/GeometricSolvers.en.md) * [Symbolic Vectors, Matrices and Arrays](https://reference.wolfram.com/language/guide/SymbolicArrays.en.md) * [Numerical Evaluation & Precision](https://reference.wolfram.com/language/guide/NumericalEvaluationAndPrecision.en.md) * [Mesh-Based Geometric Regions](https://reference.wolfram.com/language/guide/MeshRegions.en.md) * [`Polygons`](https://reference.wolfram.com/language/guide/Polygons.en.md) * [`Polyhedra`](https://reference.wolfram.com/language/guide/Polyhedra.en.md) ## Related Links * [NIntegrate Explorer](https://reference.wolfram.com/language/GUIKit/tutorial/NIntegrateExplorer.en.md) * [NKS\|Online](http://www.wolframscience.com/nks/search/?q=NIntegrate) * [A New Kind of Science](http://www.wolframscience.com/nks/) ## History * Introduced in 1988 (1.0) \| Updated in 1996 (3.0) ▪ 2003 (5.0) ▪ [2007 (6.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn60.en.md) ▪ [2010 (8.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn80.en.md) ▪ [2014 (10.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn100.en.md)