Wolfram Research, Inc.

3.9.3 Numerical Integration

Numerical integration functions.

This finds a numerical approximation to the integral .

In[1]:= NIntegrate[Exp[-x^3], {x, 0, Infinity}]

Out[1]=

Here is the numerical value of the double integral .

In[2]:= NIntegrate[x^2 + y^2, {x, -1, 1}, {y, -1, 1}]

Out[2]=

An important feature of NIntegrate is its ability to deal with functions that "blow up" at known points. NIntegrate automatically checks for such problems at the end points of the integration region.

The function blows up at , but NIntegrate still succeeds in getting the correct value for the integral.

In[3]:= NIntegrate[1/Sqrt[x], {x, 0, 1}]

Out[3]=

Mathematica can find the integral of exactly.

In[4]:= Integrate[1/Sqrt[x], {x, 0, 1}]

Out[4]=

NIntegrate detects that the singularity in at is not integrable.

In[5]:= NIntegrate[1/x, {x, 0, 1}]

Out[5]=

NIntegrate does not automatically look for singularities except at the end points of your integration region. When other singularities are present, NIntegrate may not give you the right answer for the integral. Nevertheless, in following its adaptive procedure, NIntegrate will often detect the presence of potentially singular behavior, and will warn you about it.

NIntegrate does not handle the singularity in in the middle of the integration region. However, it warns you of a possible problem. In this case, the final result is numerically quite close to the correct answer.

In[6]:= NIntegrate[1/Sqrt[Abs[x]], {x, -1, 2}]

Out[6]=

If you know that your integrand has singularities at particular points, you can explicitly tell NIntegrate to deal with them. NIntegrate[expr, x, xmin, , , ... , xmax] integrates expr from xmin to xmax, looking for possible singularities at each of the intermediate points .

This again gives the integral , but now explicitly deals with the singularity at .

In[7]:= NIntegrate[1/Sqrt[Abs[x]], {x, -1, 0, 2}]

Out[7]=

You can also use the list of intermediate points in NIntegrate to specify an integration contour to follow in the complex plane. The contour is taken to consist of a sequence of line segments, starting at xmin, going through each of the , and ending at xmax.

This integrates around a closed contour in the complex plane, going from , through the points , 1 and , then back to .

In[8]:= NIntegrate[1/x, {x, -1, -I, 1, I, -1}]

Out[8]=

The integral gives , as expected from Cauchy's Theorem.

In[9]:= N[ 2 Pi I ]

Out[9]=

There are a number of ways you can control the operation of NIntegrate. First, you may want to specify the accuracy of the answers you are trying to get. If you use Integrate, and then apply N to get numerical results, then the second argument you give to N specifies the precision to use in the internal computations used to do the numerical integral. You can also specify this precision using the option WorkingPrecision -> n for NIntegrate.

You should realize, however, that the option WorkingPrecision specifies only the precision to use for internal computations done by NIntegrate. The final answer that NIntegrate gives will almost always have a lower precision. You can nevertheless use NIntegrate to try and get an answer with a particular precision or accuracy by setting the options PrecisionGoal or AccuracyGoal. The default setting PrecisionGoal -> Automatic attempts to get answers with a precision equal to WorkingPrecision minus 10 digits. In general, NIntegrate continues until it achieves either of the goals specified by PrecisionGoal or AccuracyGoal.

This evaluates using 40-digit precision for internal computations. The precision goal in this case is automatically set to 30 digits.

In[10]:= NIntegrate[1/x, {x, 1, 2}, WorkingPrecision->40]

Out[10]=

The difference from the exact result is consistent with 0.

In[11]:= % - N[Log[2], 40]

Out[11]=

In fact, the result was correct to about 33 digits.

In[12]:= Accuracy[%]

Out[12]=

Options for NIntegrate.

When NIntegrate tries to evaluate a numerical integral, it samples the integrand at a sequence of points. If it finds that the integrand changes rapidly in a particular region, then it recursively takes more sample points in that region. The parameters MinRecursion and MaxRecursion specify the minimum and maximum number of levels of recursive subdivision to use. Increasing the value of MinRecursion guarantees that NIntegrate will use a larger number of sample points. MaxRecursion limits the number of sample points which NIntegrate will ever try to use. Increasing MinRecursion or MaxRecursion will make NIntegrate work more slowly. SingularityDepth specifies how many levels of recursive subdivision NIntegrate should try before it concludes that the integrand is "blowing up" at one of the endpoints, and does a change of variables.

With the default settings for all options, NIntegrate misses the peak in near , and gives the wrong answer for the integral.

In[13]:= NIntegrate[Exp[-x^2], {x, -1000, 1000}]

Out[13]=

With the option MinRecursion->3, NIntegrate samples enough points that it notices the peak around . With the default setting of MaxRecursion, however, NIntegrate cannot use enough sample points to be able to expect an accurate answer.

In[14]:= NIntegrate[Exp[-x^2], {x, -1000, 1000},

MinRecursion->3]

Out[14]=

With this setting of MaxRecursion, NIntegrate can get an accurate answer for the integral.

In[15]:= NIntegrate[Exp[-x^2], {x, -1000, 1000},

MinRecursion->3, MaxRecursion->10]

Out[15]=

Another way to solve the problem is to make NIntegrate break the integration region into several pieces, with a small piece that explicitly covers the neighborhood of the peak.

In[16]:= NIntegrate[Exp[-x^2], {x, -1000, -10, 10, 1000}]

Out[16]=

For integrals in many dimensions, it can take a long time for NIntegrate to get a precise answer. However, by setting the option MaxPoints, you can tell NIntegrate to give you just a rough estimate, sampling the integrand only a limited number of times.

This gives an estimate of the volume of the unit sphere in three dimensions.

In[17]:= NIntegrate[If[x^2 + y^2 + z^2 < 1, 1, 0], {x, -1, 1},

{y, -1, 1}, {z, -1, 1}, MaxPoints->10000]

Out[17]=

Here is the precise result.

In[18]:= N[4/3 Pi]

Out[18]=