Integrals

Indefinite integrals

Indefinite integral of a function is a primitive (original, antiderivative) of this function. If

(∂_x)F[x]=f[x],

then the indefinite integral of f[x] is

F[x]=∫f[x]dx,

Indefinite integral is just the inverse to derivative. It is defined up to an additive constant. Indefinite integrals of some elementary functions can be analytically expressed via elementary functions. In simple cases it is sufficient to know a derivarive to obtain an integral

Power function:

(∂_x)(x^α)

∫(x^α)\[DifferentialD]x

Note that Mathematica does not output an additive constant in indefinite integrals. Check:

(∂_x)((x^(1+α))/(1+α))

OK.

Trigonometric functions:

∫Cos[x]\[DifferentialD]x

Check:

(∂_x)Sin[x]

OK.

Inverse trigonometric functions:

∫(1/(1+(x^2)))\[DifferentialD]x

(∂_x)ArcTan[x]

In some cases analytical calculation of indefinite integrals is more complicated and requires more work and knowledge of tricks.

∫(1/(Sin[x]))\[DifferentialD]x

This is

Log[Tan[(x/2)]]

Check

Simplify[(∂_x)Log[Tan[(x/2)]]]

but it is not obvious how to force Mathematica to bring the result in this form. The idea of the calculation must be using Tan[(x/2)] as the integration variable.

y= Tan[(x/2)],      dy=((∂_x)Tan[(x/2)])dx =(1/2) ((Sec[(x/2)])^2)dx

(∂_x)Tan[(x/2)]

Now

∫(1/(Sin[x]))dx=∫((2((Cos[(x/2)])^2))/(Sin[x]))dy=∫((2((Cos[(x/2)])^2))/(2Sin[(x/2)]Cos[(x/2)]))dy=∫Cot[(x/2)]dy=∫(1/y)dy=Log[y]=Log[Tan[(x/2)]].

Another example

∫(1/(x(1+(x^2))))\[DifferentialD]x

This integral can be calculated by splitting the fraction in two:

Apart[(1/(x(1+(x^2))))]

and then using y=(x^2) as the integration variable in the second term. Probably it does not make a lot of sense to spend much time on analytical calculating integrals by hand because it can be very tedious and computer does it fast.

∫((1-2x+3(x^2))/(x(1+(x^2))(1+x+2(x^2))))\[DifferentialD]x

Although derivatives of elementary functions are elementary functions, primitives of elementary functions can be special functions, or analytical results in terms of known functions do not exist. Below are so-called elliptic integrals of two different kinds

∫(1/(Sqrt[1-m ((Sin[x])^2)]))\[DifferentialD]x
∫(Sqrt[1-m ((Sin[x])^2)])\[DifferentialD]x

In particular, EllipticE[x, 1] can be reduced to sinusoidal functions

Plot[{EllipticE[x,1],Sin[x]},{x,0,2π},PlotStyle->{Black,{Black,Dashed}}]

but, in general, elliptic integrals are special functions, non-reducible to elementary functions. Another example of a special function is the error function Erf[x]:

∫(E^(-(x^2)))\[DifferentialD]x

Plot[Erf[x],{x,-5,5}]

In other cases indefinite integrals cannot be expressed via any known function and have to be calculated numerically.

∫((E^(-(x^2)))/(1+(x^2)))\[DifferentialD]x

Definite integrals

Definition

Definite integrals are defined as the limit of the sum

(∫_a^b)f[x]dx=(lim_(N->∞))(∑_(n=0)^(N-1))f[(Subscript[x, n])]Δx,       (Subscript[x, n])=a+(b-a)(n/N),    Δx≡h=((b-a)/N),

as illustrated in the figures below.

f[x_]:=1+(x^2);
NN=5;
a=0;   b=1;
xn[n_]:=a+(b-a)(n/NN)

Plotting with lines

Pf=Plot[f[x],{x,a,b},PlotRange->{0,2},PlotStyle->{Thick}];
LinesVert=Graphics[Table[Line[{{xn[n],0},{xn[n],Max[f[xn[n]],f[xn[n-1]]]}}],{n,0,NN}]];
LinesHoriz=Graphics[Table[Line[{{xn[n],f[xn[n]]},{xn[n+1],f[xn[n]]}}],{n,0,NN-1}]];
Show[Pf,LinesVert,LinesHoriz]

Plotting with rectangles

Pf=Plot[f[x],{x,a,b},PlotRange->{0,2},PlotStyle->{Thick,Black}];
Rectangles=Graphics[Table[{EdgeForm[{Thick, Dashed}],Pink,Rectangle[{xn[n],0},{xn[n+1],f[xn[n]]}]},{n,0,NN-1}]];
Show[Pf,Rectangles,Pf]

Geometrical interpretation of a difinite integral is the area below the curve representing the function f[x]. The exact area is first approximated by the total area of the rectangles and then the limit if infinitesimaly small rectangles is taken.

For increasing as above, the sum always underestimates the definite integral, whereas for decreasing functions it overestimates the definite integral.

f[x_]:=1-(x^2);
NN=5;
a=0;   b=1;
xn[n_]:=a+(b-a)(n/NN)
Pf=Plot[f[x],{x,a,b},PlotRange->{0,1},PlotStyle->{Thick,Black}];
Rectangles=Graphics[Table[{EdgeForm[{Thick, Dashed}],Pink,Rectangle[{xn[n],0},{xn[n+1],f[xn[n]]}]},{n,0,NN-1}]];
Show[Pf,Rectangles,Pf]

Note that there is no left-right symmetry in this definition of the integral and the rectangles can be chosen in two different ways, over- and underestimating the integral.

f[x_]:=1-(x^2);
NN=10;
a=0;   b=1;
xn[n_]:=a+(b-a)(n/NN)
Pf=Plot[f[x],{x,a,b},PlotRange->{0,1},PlotStyle->{Thick,Black}];
Rectangles1=Graphics[Table[{EdgeForm[{Thick, Dashed}],Pink,Rectangle[{xn[n],0},{xn[n+1],f[xn[n]]}]},{n,0,NN-1}]];
Rectangles2=Graphics[Table[{EdgeForm[{Thick, Dashed}],LightBlue,Rectangle[{xn[n],0},{xn[n+1],f[xn[n+1]]}]},{n,0,NN-1}]];
Show[Pf,Rectangles1,Rectangles2,Pf]

Nevertheless, the limit N->∞ in both methods is the same.

Numerical calculation of integrals

Quadrature formulas

Let us discuss the convergence of the sums (Subscript[I, N]) to the integral I. The difference (Subscript[ΔI, N])=|(Subscript[I, N])-I| is due to the triangles located below the curve above the rectangles. The area of these triangles can be estimated as

(Subscript[ΔI, N])=(1/2)(∑_(n=0)^(N-1))f'[(Subscript[x, n])](((Δx))^2)~(((Δx))^2)N~Δx=h

This is the so-called first-order convergence that is slow and insatisfactory for numerical approximation of an integral by a sum. Convergence can be strongly improved if one replaces rectangles by trapezias.

f[x_]:=Cos[x];
NN=4;
a=0;   b=π/2;
xn[n_]:=a+(b-a)(n/NN)
Pf=Plot[f[x],{x,a,b},PlotRange->{0,1},PlotStyle->{Thick,Black}];
Trapezias=Graphics[Table[{EdgeForm[{Thick, Dashed}],Pink,Polygon[{{xn[n],0},{xn[n],f[xn[n]]},{xn[n+1],f[xn[n+1]]},{xn[n+1],0}}]},{n,0,NN-1}]];
Show[Pf,Trapezias]

This construction has the right-left symmetry, and the approximation of the integral has the form

(∫_a^b)f[x]dx=(lim_(N->∞))(∑_(n=0)^(N-1))(1/2)(f[(Subscript[x, n])]+f[(Subscript[x, n+1])])Δx,       (Subscript[x, n])=a+(b-a)(n/N),    Δx≡h=((b-a)/N),

or

(∫_a^b)f[x]dx=((1/2)f[(Subscript[x, 0])]+(∑_(n=1)^(N-1))f[(Subscript[x, n])]+(1/2)f[(Subscript[x, N])])Δx +O((h^2)).

This approximation of the integral by the sum, the so-called trapezoidal rule, is one of many quadrature formulas for integrals. The difference of the trapezoidal rule and the rectangle formula used for the definition of the integral is only in the coefficients in front of the boundary terms, here 1/2. It can be shown that the convergence of the trapezoidal frule is (h^2) that is much better that the convergence of the rectangle formula.

Generally, an integral can be represented by a sum of the form

(∫_a^b)f[x]dx=(∑_(n=0)^N)(Subscript[w, n])f[(Subscript[x, n])]

where (Subscript[x, n]) are nodes and (Subscript[w, n]) are weights. Optimizing positions of the nodes and values of the weights, very precise and fast convergent quandatures can be constructed. In all cases the function is approximated by polynomials within the interval between nodes and these polynomials are integrated resulting in a quadrature formula. One group of quadratures are those with equidistant nodes, such as the rectangle and trapezoidal rules above. Approximating the function by quadratic or cubic polynomials, one obtains more precise Simpson formulas in which weights are alternating. Durand's rule

(∫_a^b)f[x]dx=((2/5)f[(Subscript[x, 0])]+(11/10)f[(Subscript[x, 1])]+(∑_(n=2)^(N-2))f[(Subscript[x, n])]+(11/10)f[(Subscript[x, N])]+(2/5)f[(Subscript[x, N])])Δx .

is only slightly more complicated than the trapezoidal rule, all the weight except those at the four boundary points being equal to Δx. Durand's rule is much more precise and fast converging than the trapezoidal rule. The order of its convergence is difficult to find but one can investigate it practically comparing convergence of the rectangular, trapezoidal, and Durand rules.

f[x_]:=Cos[x];
a=0;   b=π/2.;
II=(∫_a^b)f[x]\[DifferentialD]x    
xNn[NN_,n_]:=a+(b-a)(n/NN)
NMax=50;

INRectangular[NN_]:=((b-a)/NN)(∑_(n=1)^NN)f[xNn[NN,n]]

INTrapezoidal[NN_]:=((b-a)/NN)((1/2)f[xNn[NN,0]]+(∑_(n=1)^(NN-1))f[xNn[NN,n]]+(1/2)f[xNn[NN,NN]])
INDurand[NN_]:=((b-a)/NN)((2/5)f[xNn[NN,0]]+(11/10)f[xNn[NN,1]]+(∑_(n=2)^(NN-2))f[xNn[NN,n]]+(11/10)f[xNn[NN,NN-1]]+(2/5)f[xNn[NN,NN]])

INRectangularList=Table[{NN,INRectangular[NN]},{NN,2,NMax}];
INTrapezoidalList=Table[{NN,INTrapezoidal[NN]},{NN,2,NMax}];
INDurandList=Table[{NN,INDurand[NN]},{NN,2,NMax}];

P1=Plot[II,{NN,3,1.03NMax},PlotRange->{0.97,1.001},PlotStyle->{Black,Thick},AxesLabel->{"N","(Subscript[I, N])"},BaseStyle->{FontSize->14,FontFamily->"Times",Plain}];
PINRectangularList=ListPlot[INRectangularList,PlotStyle->{Green,PointSize[0.015]}];
PINTrapezoidalList=ListPlot[INTrapezoidalList,PlotStyle->{Blue,PointSize[0.015]}];
PINDurandList=ListPlot[INDurandList,PlotStyle->{Red,PointSize[0.015]}];

Show[P1,PINRectangularList,PINDurandList,PINTrapezoidalList]

To better illustrate convergence of different quadrature rules, one can make a Log-Log plot of (Subscript[ΔI, N])=|(Subscript[I, N])-I|. Our expectation is

(Subscript[ΔI, N])=a (h^η)=b (N^(-η)),

where η is the order of convergence (η = 1 for the rectangular rule). In the Log-Log form this relation is a straight line:

Log[(Subscript[ΔI, N])]=Log[b] -η Log[N],

and η can be found from its slope.

NMax=100;
ΔINRectangularList=Table[{NN,II-INRectangular[NN]},{NN,2,NMax}];
ΔINTrapezoidalList=Table[{NN,II-INTrapezoidal[NN]},{NN,2,NMax}];
ΔINDurandList=Table[{NN,Abs[II-INDurand[NN]]},{NN,2,NMax}];
ListLogLogPlot[{ΔINRectangularList,ΔINTrapezoidalList,ΔINDurandList},PlotRange->All,PlotStyle->{Green,Blue,Red},AxesLabel->{"N","(Subscript[ΔI, N])"},BaseStyle->{FontSize->14,FontFamily->"Times",Plain}]

Indeed, all three dependences are asymptotically straight lines. Surprizingly, the slopes for the trapezoidal and Durand rules are the same, that is, their order is η = 2. However, the coefficient b is by a factor 8 or so smaller for the Durand rule.

There are even more precise quadratures with equal weights for most of the nodes, except for the boundary terms. If the function is exponentially small near the ends of the interval, contributions of the boundary terms become negligible and one can use the simplest quadrature with all equal weights to achieve exponentially precise approximations for the integral. As an example we can calculate the integral

(1/(Sqrt[π]))(∫_(-∞)^∞)(E^(-(x^2)))\[DifferentialD]x

Plot[(1./(Sqrt[π]))(E^(-(x^2))),{x,-6,6}]

numerically using the simplest quadrature rule

(∫_a^b)f[x]dx=(∑_(n=0)^N)f[(Subscript[x, n])]

f[x_]:=(1./(Sqrt[π]))(E^(-(x^2)));
a=-6;   b=6;  (* outside this interval the function becomes very small *)

xNn[NN_,n_]:=a+(b-a)(n/NN)
NMax=20;

IN[NN_]:=((b-a)/NN)(∑_(n=0)^NN)f[xNn[NN,n]]
INList=Table[{NN,IN[NN]},{NN,2,NMax}]
P1=Plot[1,{NN,3,NMax},PlotRange->{0.5,1.5},PlotStyle->{Black,Thick},AxesLabel->{"N","(Subscript[I, N])"},BaseStyle->{FontSize->14,FontFamily->"Times",Plain}];
PINList=ListPlot[INList,PlotStyle->{Red,PointSize[0.015]}];
Show[P1,PINList]

Here the convergence is faster than everything we have seen before. To better see it, let us make a Log-Log plot of (Subscript[ΔI, N])

ΔINList=Table[{NN,Abs[II-IN[NN]]},{NN,2,NMax}];
ListLogLogPlot[ΔINList,AxesLabel->{"N","(Subscript[ΔI, N])"},BaseStyle->{FontSize->14,FontFamily->"Times",Plain}]

This is not a straight line and the values of (Subscript[ΔI, N]) become extremely small. The results suggest an exponential dependence of (Subscript[ΔI, N]) on N. We will see that this dependence has the form

(Subscript[ΔI, N]) =b Exp[-c (N^2)]

so that

Log[(Subscript[ΔI, N])] =Log[b]-c (N^2).

Let us thus plot Log[(Subscript[ΔI, N])] vs (N^2) that is expected to be a straight line.

ΔINList2=Table[{(NN^2),Abs[II-IN[NN]]},{NN,2,NMax}];
ListLogPlot[ΔINList2,AxesLabel->{"(N^2)","(Subscript[ΔI, N])"},BaseStyle->{FontSize->14,FontFamily->"Times",Plain}]

Returning to more sophisticated quadratures, one can choose positions of the nodes in a way to achieve a faster convergence. The corresponding quadratures with non-equidistant nodes have been proposed by Gauss and they are superior to the quadratures with equidistant nodes.

Since highly efficient quadratures are implemented in the software, it does not make sense to study them in detail in our course.

Numerical integration with Mathematica

Mathematica computes integrals numerically with the help of NIntegrate command

NIntegrate[((E^(-(x^2)))/(1+(x^2))),{x,-∞,∞}]

that is equivalent to

N[(∫_(-∞)^∞)((E^(-(x^2)))/(1+(x^2)))\[DifferentialD]x]

This particular integral can be calculated analytically via an error function

(∫_(-∞)^∞)((E^(-(x^2)))/(1+(x^2)))\[DifferentialD]x
N[%]

but if the limits are finite, numerical evaluation is the only option. For numerical calculation, Mathematica can use a lot of different quadratures. Using NIntegrate, one can specify which method should be used as an option. At least for smooth functions as above, there is no big difference between the methods.

Problems in numerical integration can arise if the integrand is singular

NIntegrate[(1/((Abs[(x-1)(x-1.5)(x-2.3)(x-3)])^(3/4))),{x,0,∞}]

Here the imaginary part is an artefact and the integrator complains about recursions. (Integral is computed recursively adapting the number of integration points) After adding the option for maximal number of recursions it becomes better

NIntegrate[(1/((Abs[(x-1)(x-1.5)(x-2.3)(x-3)])^(3/4))),{x,0,∞},MaxRecursion->100]

although the imaginary part is still there. Specifying positions of the singularities removes the imaginary part

NIntegrate[(1/((Abs[(x-1)(x-1.5)(x-2.3)(x-3)])^(3/4))),{x,0,1,1.5,2.3,3,∞}]

NIntegrate[(1/((Abs[(x-1)(x-1.5)(x-2.3)(x-3)])^(3/4))),{x,0,1,1.5,2.3,3,∞},MaxRecursion->100]

In this case the result was not bad in all cases. But the problems with singularities can become more severe in the case of multi-dimentional integrals.

Another example is integration of a strongly oscillating function Sin[(1/(x^2))]

Plot[Sin[(1/(x^2))],{x,0,10},PlotRange->All]

Straightforward numerical integration yields

NIntegrate[Sin[(1/(x^2))],{x,0,∞}]
NIntegrate[Sin[(1/(x^2))],{x,0,∞},MaxRecursion->100]

The accuracy in both cases is insufficient, as follows from Mathematica warnings and comparison with the exact value of this integral

Integrate[Sin[(1/(x^2))],{x,0,∞}]
N[%]

In such cases is may be useful to analytically transform the integral before numerical integration. With the new variable y=1/(x^2) one obtains

x=(1/(Sqrt[y])),    dx=-(1/2)(dy/(y^(3/2))),   I=(∫_0^∞)Sin[(1/(x^2))]dx=(1/2)(∫_0^∞)((Sin[y])/(y^(3/2)))dy

and

NIntegrate[(1/2)((Sin[y])/(y^(3/2))),{y,0,∞}]

without any numerical problems. Indeed, the integrand ((Sin[y])/(y^(3/2))) is less oscillatory than the original Sin[(1/(x^2))]:

Plot[(1/2)((Sin[y])/(y^(3/2))),{y,0,30},PlotRange->{-0.1,0.5}]

Analytical calculation of definite integrals

If the indefinite integral F[x] of f[x] is known, the definite integral of f[x] is given by the well-known formula

(∫_a^b)f[x]dx=F[b]-F[a].

If F[x] cannot be obtained, still in some cases definite integrals, especially with infinite limits, can be analyticaly calculated by special methods or tricks. For instance, the integral arising in statistical thermodynamics

(∫_0^∞)(x/((E^x)-1))\[DifferentialD]x

can be reduced to the series

(∑_(n=1)^∞)(1/(n^2))

by expanding the integrand into the geometrical progression and integrating the resulting terms as

(∫_0^∞)x (E^(-n x))\[DifferentialD]x

One obtains

(∫_0^∞)(x/((e^x)-1))dx=(∫_0^∞)((x (e^(-x)))/(1-(e^(-x))))dx=(∫_0^∞)x((e^(-x))+(e^(-2x))+(e^(-3x))+...)dx=1+(1/(2^2))+(1/(3^2))+...=((π^2)/6)

Summation of this series is more involved and it can be done by expressing it via a double integral that can be analytically calculated.

The most powerful method of calculating definite integrals is using complex calculus.

If the integrand depends on parameters, the possibility to calculate the definite integral depends of these parameters. Thus the result includes conditions on the parameters as in the integral above. To obtain a nicer output in the cases when the convergence conditions are clear, one can use the option Assumptions in Integrate:

Integrate[x (E^(-n x)),{x,0,∞},Assumptions->n>0]

In this way one can define a function like

fn[n_]:=Integrate[x (E^(-n x)),{x,0,∞},Assumptions->n>0]
fn[n]

that is impossible without using Assumptions

Convergence of integrals

Integrals can be convergent or divergent at singular points and at infinity. The most frequent type of singularity is a power singularity (((x-a))^α). The integral can also contain a non-singular function having a finite value at the singularity. Integral

(∫_a^b)(((x-a))^α)f[x]dx

converges if α > -1 independently of the exact form of f[x] that can be replaced by f[a] for the investigation of convergence. For -1 < α < 0 the integrand diverges and numerical calculation becomes slightly more difficult. To avoid numerical calculating a diverging function, one can rewrite the integral in the form

(Subscript[(∫_a^b)(((x-a))^α)f[x]dx=f[a](∫_a^b)(((x-a))^α)+(∫_a^b)(((x-a))^α)(f[x]-f[a])dx=((f[a])/(α+1))(((x-a))^α)|, a]^b)+(∫_a^b)(((x-a))^α)(f[x]-f[a])dx=((f[a])/(α+1))(((b-a))^α)+(∫_a^b)(((x-a))^α)(f[x]-f[a])dx.

Here the singular part of the integral has been calculated analytically and the remainder is an integral of a non-diverging function. Indeed, near x=a one has

(((x-a))^α)(f[x]-f[a])≈(((x-a))^α)f'[a](x-a)~(((x-a))^(α+1)).

This method can be used in many cases to kill singularities. If the integrand behaves as a power at infinity,

(∫_a^∞)(x^α)f[x]dx,     f[∞]=const,

the integral converges if α < -1, independently of the exact form of f[x]. By comparison with this type of integrals, one can investigate convergence of series.