Calculus

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

In this notebook we present some of the basic techniques for performing standard operations in calculus, including differentiation, Taylor expansion, and integration.  Both symbolic and numerical integration are covered and brief discussions of line and contour integration are given also.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Needs["Utilities`Notation`"]

Differentiation

Partial derivatives

Partial derivatives are produced by the function D[f,x], where the expression f should be either an implicit or explicit function of the variable x for which the derivative is sought.  Partial derivatives with respect to several variables simultaneously are obtained by using a sequence of variables in the argument list.  Multiple derivatives with respect to the same variable are specified by arguments of the form {x,n} where n is the degree of differentiation or by including n instances of x in the sequence of differentiation variables.

There are many equivalent input forms for this function to reflect common usages.  The prime notation is convenient for writing ordinary differential equations and can be entered directly from the keyboard; a slightly more appealing form is obtained by entering Esc'Esc in a superscript box.  The (∂_({x, n}, {y, m}))f[x, y] notation is more flexible and more rigorous, superseding the traditional \(TraditionalForm\`∂\^\(n + m\)f[x, y]\/\(∂\^n x ∂\^m y\)\) notation with its cumbersome fractions which do not divide out.  That traditional notation cannot be used without constructing special notation palettes first.  The partial derivative symbol is entered as EscpdEsc.

D[f[x], 
        x] = ∂f[x]/∂x = f'[
            x] = (f^′)[
              x] = (∂_x) f[x]

\(D[f[x], {x, n}] = 
      \(∂\^n f[x]\/∂x\^n = ∂\_{x, n}f[x]\)\)

D[f[x, y], x, y] = 
      (∂^2) f[x, y]/(∂x ∂y) = 
        (∂_(x, y))f[x, y]

D[f[x, y], {x, n}, {y, m}] = 
      (∂^(n + m))f[x, y]/((∂^n) x 
            (∂^m) y) = (∂_({x, n}, {y, m}))f[x, y]

D[f[x], x]

D[f[x, y], {x, n}, {y, m}]

(∂_({x, n}, {y, m}))f[x, y]

(∂_(x, x))f[x]

\(∂\_{x, 4}f[x]\)

f''''[x]

Note that entering

(f^((4)))[X]

with an exponent box returns a power instead of a derivative, parentheses notwithstanding.  You would need to develop special notation palettes to use traditional notations like (f^((4)))[x] or (∂^4) f/∂(x^4).

D employs all the basic rules for differentiating products and quotients and uses the chain rule.

Clear[f, g]

D[f[x] g[x], x]

D[f[x y], x]

D[f[g[x]], x]

D[(∫_a^x) f[t] \[DifferentialD]t, x]

D[(∫_a^(Cos[x]))f[t] \[DifferentialD]t, x]

It is always important to remember that because D evaluates partial derivatives it treats all symbols that are not declared variable as constants.  Symbols are declared as variable either by inclusion in the argument list or by means of the option NonConstants->list.  Suppose that x is an implicit function of t, where x is the position of a parcel of fluid and t is time.  The expression

D[f[x, t], t]

then evaluates the derivative of f with respect to t at a fixed position, whereas

D[f[x, t], t, NonConstants -> x]

evaluates the convective derivative for a particular moving parcel of fluid.  Of course, the same effect is achieved by declaring x to be an explicit function of t.

D[f[x[t], t], t]

Locate the extrema of (x^3) - 2  (x^2) - 10  x + 5.

Find the maximum of ((x^4)/((e^x) + 1)) for positive x (both position and value).  [Hint: simplify the equation first to eliminate the trivial but spurious solution.  Then use Solve, NSolve, or FindRoot; all should agree!]

Jacobian for coordinate transformations

J[u, v, w]=(({{∂(f/∂)u, ∂(f/∂)v, ∂(f/∂)w}, {∂(g/∂)u, ∂(g/∂)v, ∂(g/∂)w}, {∂(h/∂)u, ∂(h/∂)v, ∂(h/∂)w}}))

Total derivatives

The total derivative Dt[f,x] evaluates the derivative ((d,f)/(d,x)) by interpreting any variable that appears in the expression f as a function of x.  Consider, for example, a function of position x and time t where x is an implicit function of time.  The total derivative with respect to time

Dt[f[x, t], t]

contains two contributions, one arising from the explicit dependence upon t and a second arising from the implicit dependence of x upon t.  Similarly, multiple or mixed derivatives can be defined also.

Dt[f[x, t], x, x, t]

A closely related quantity is the total differential

Dt[f[x, y, z, t]]

which represents the differential change in a function arising from differential changes in each variable that is not explicitly constant.  One or more variables can be held constant by using the option Constants to specify a list of constants.

Dt[f[x, y, z, t], Constants -> {x, y}]

Alternatively, we can eliminate selected differentials

Dt[f[x, y, z, t]] /. {Dt[x] -> 0, Dt[y] -> 0}

and obtain a slightly simpler result that does not propagate the option in the result.  The best choice depends upon the application.

However, if expressions involve symbols which should always be constant, such as constants of nature, it is best to declare them constants by setting the appropriate attribute.

SetAttributes[{ℏ, m}, Constant]

\(Dt[\(-\(ℏ\^2\/\(2  m\)\)\) ∂\_{x, 2}ψ[x, t], 
      t]\)

Example: particle in uniform gravity

The energy of a particle in uniform gravity has the form

ClearAll[m, g, x, y, z, t];
energy  =  
      m g z[t]  +  (1/2) 
          m ((x^′)[t]^2 + (y^′)[t]^2 + (z^′)[t]^2)

The time dependence of the energy is given by the total derivative

Dt[energy, t]

However, this result includes derivatives of mass with respect to time which should not appear for a simple particle that neither sheds nor accretes mass.  We specify that m and g are constants using the option Constants, such that

Dt[energy, t, Constants -> {m, g}]

Furthermore, Newton's second law can be used to simplify this result

Dt[energy, t, 
        Constants -> {m, g}] /. {(x^″)[t] -> 
          0, (y^″)[t] -> 
          0, (z^″)[t] -> -g}

The dependence of the energy on height is then seen to vanish; in other words, energy is conserved.

Fluid flows at rate \(TraditionalForm\`Q = \(π\ p\ r\^4\)\/\(8  ℓ\ η\)\) through a cylindrical tube with radius r and length ℓ subject to a pressure differential p.  The coefficient of viscosity is η.  Suppose that the length of the tube is constant but that the radius expands at rate \(TraditionalForm\`r\& . \) as the pressure differential decreases at rate \(TraditionalForm\`p\& . \).  Find an expression for the time rate of change for Q.  Under what conditions is Q constant?

Implicit differentiation

Suppose that two variables, x and y, are connected by an equation of constraint, eq.  Evaluation of Dt[eq] will then produce a relationship between variations of x, Dt[x], and variations of y, Dt[y].  For example, suppose that a ladder of length ℓ leans against a wall but that the base of the ladder slides along the ground.  To compare the rate the top falls with the rate the bottom slides, we differentiate the equation

eq  =  (x^2) + (y^2) == (ℓ^2);
newEq = Dt[eq]

that evaluates the length of the ladder in terms of the coordinates of its ends.  Requiring the length to be constant then implies

sol = Solve[newEq /. Dt[ℓ] -> 0, 
          Dt[y]][[1]]

We can express the result in terms of the position of the bottom ladder by using the solution to the constraint equation

sol[[1, 
        2]]  /. Solve[eq, 
          y][[2]]

Although we could have solved the equation of constraint first and then differentiated the result, we chose the approach above in order to illustrate the technique of implicit differentiation with a simple example.

A lemniscate is described by the equation ((((x^2) + (y^2)))^2) == (x^2) - (y^2).  Use ImplicitPlot to display this curve and implicit differentiation to find the points at which its tangent is horizontal.

A person who is 6 ft tall walks toward a street lamp that is 20 ft above its base.  How fast does that person's shadow shrink?  Use similar triangles to establish a constraint equation and implicit differentiation to relate the length change to the distance walked.

Derivatives as functionals

D will evaluate the derivatives for any function for which Mathematica knows the necessary rules, but when rules are not available it returns a generic Derivative object.

g = (∂_({x, n}, {y, m}))f[x, y]

FullForm[g]

Notice that Derivative[n,m][f] represents a pure function whose arguments are [x,y].

Head[g] // FullForm

(f^′) // FullForm

This internal representation of derivatives actually embodies a very important concept.  Consider for a moment the traditional notation f'[x].  This notation really means

(f^′)[
        x] =  (Subscript[((∂f[t]/∂t)), t = x])

where t is a dummy variable that replaced by the desired argument after the derivative is evaluated.  Hence, we interpret (f^′)[2] as (Subscript[(((f^′)[x])), x = 2]) even though no variable appears in the traditional notation.  Thus, differentiation is an operator which acts upon a function to produce another function.  In Mathematica the result of differentiation is represented by a pure function that can be evaluated either for symbolic or numerical arguments.  The operation of differentiation itself can be viewed as a pure function whose argument is a function.

Head[g]

Head[Head[g]]

By this means one can view the symbolic differentiation rules associated with standard mathematical functions.

f = (Sin^′)

The result is a pure function which can be applied to any argument.

{f[x], f[(π/2)]}

Integrate[f[x], x]

Locating minima

Using FindMinimum

The extrema of a function (maxima, minima, and inflection points) can be determined by solving equations, either symbolically or numerically, that set first derivatives to zero and then ascertaining the characteristics of the second derivatives.  Symbolic methods will usually reveal all extrema, while numerical methods search for extrema near a specified starting point.  Mathematica provides a built-in function FindMinimum which searches for a local minimum near the starting point automatically, without the user having to establish the equations for the first derivative or to test the second.  The result is usually the nearest minimum, which may not be the absolute minimum of the function — it is simply a nearby local minimum.  If you need an absolute minimum, it is usually necessary to examine the function graphically in order to supply a reasonably good starting value or set of starting values.  You also need to check that the result truly is the desired minimum — sometimes FindMinimum wanders away from your starting values and sometimes it returns incorrect results.

Example: a function of one variable

The basic syntax is FindMinimum[f,{x, (Subscript[x, 0])}] where f  is a function of x and (Subscript[x, 0]) is the starting value.  The function must evaluate to a real number given a numerical value for x.  You can limit the range of x by including xmin, xmax after (Subscript[x, 0]).

Clear[f];
f = (1 + x - (x^2)) Exp[-(x^2)/4];

Plot[f, {x, -5, 5}];

The positive minimum is found using an appropriate choice of x.  The result is returned as a list whose first element is the value of the function as the minimum and the second is a list of replacement rules specifying the location of the minimum.

FindMinimum[f, {x, 3}]

Both minima can be found by mapping FindMinimum onto a list of starting points.

FindMinimum[f, {x, #}] &  /@ {-3, 3}

How can you locate the maximum using FindMinimum?  Try it!

Radial maxima for hydrogenic wave functions

Example: a function of two variables

The basic syntax is FindMinimum[f,{x, (Subscript[x, 0])},{y, (Subscript[y, 0])}] where f  is a function of x and y and where the starting point has coordinates {(Subscript[x, 0]), (Subscript[y, 0])}.  Additional variables can be appended as necessary.  You can limit the variable ranges in the same manner as before, also.

Clear[f];
f = (x^2) Sin[x y] Exp[-(x^2)];

ContourPlot[f, {x, -2, 2}, {y, -2, 2}, 
        ContourShading -> False];

Although ContourPlot does not reveal the sign of the function directly, inspection of the definition suggests that minima are found in quadrants 2 and 4, while maxima live in quadrants 1 and 3.  Plot3D would show this, but is unnecessary.  Therefore, we give FindMinimum starting values in quadrant 2 and it reports the following minimum.

FindMinimum[f, {x, -1}, {y, 1}]

It is useful to verify this extremum by examining the function in its immediate vicinity.

ContourPlot[f, {x, -1.1, -0.9}, {y, 1.4, 1.8}, 
        Contours -> 50, ContourShading -> False];

Classifying critical points

The Taylor series expansion of a multidimensional function can be expressed in the form

f[x] ≈ 
      f[(Subscript[x, 0])] + 
        g[(Subscript[x, 0])] . (x - (Subscript[x, 0]))  +  (1/2) (x - (Subscript[x, 0])) . 
            h[(Subscript[x, 0])] . (x - (Subscript[x, 0]))  +  ⋯

where

g[(Subscript[x, 0])]  =  
      Table[∂f[x]/∂(Subscript[x, i]) 
          (Subscript[\[VerticalSeparator], x -> (Subscript[x, 0])]), {i, 1, n}]

is the gradient vector evaluated at (Subscript[x, 0]) and

h[(Subscript[x, 0])]  =  
      Table[(∂^2) f[x]/(∂(Subscript[x, i]) ∂(Subscript[x, j])) 
              (Subscript[\[VerticalSeparator] , x -> (Subscript[x, 0])]), {i, 1, n}, {j, 
          1, n}]

is the matrix of second derivatives sometimes known as the Hessian matrix.  We assume that the function f[x] is analytic at (Subscript[x, 0]) and, hence, that the Hessian matrix is symmetric.  Therefore, an important theorem of linear algebra tells us that an n-dimensional Hessian matrix has n real eigenvalues and that the corresponding eigenvectors are mutually orthogonal.  Furthermore, the Hessian matrix is diagonalized by the transformation d  =  (S^(-1)) . h . S where S is the orthogonal matrix whose columns contain the normalized eigenvectors of  h and where d is the diagonal matrix composed of the corresponding eigenvalues.

A point (Subscript[x, 0]) where the gradient vanishes, such that g[(Subscript[x, 0])] == 0, is called a critical point and is considered an ordinary critical point if it lies within the domain of analyticity of f[x].  In one dimension an ordinary critical point is a local minimum if the second derivative is positive, a local maximum if the second derivative is negative, or an inflection point if the second derivative vanishes.  In two dimensions, a critical point is classified as a saddle point if it is a local minimum for variations along one direction (not necessarily a coordinate direction) while simultaneously a local maximum for variations along the orthogonal direction.  Finally, the function exhibits a ridge if there exists a direction in which the both the first and second derivatives vanish.  Simple examples of a saddle and a ridge are shown below.

Plot3D[(y^2) - (x^2), {x, -2, 2}, {y, -2, 2}, 
        PlotLabel -> "A Simple Saddle"];

Plot3D[(y^2) (1 - (x^2)), {x, -2, 2}, {y, -2, 2}, 
        PlotLabel -> "A Simple Ridge", 
        ViewPoint -> {-0.122,  -3.114,  1.260}];

More generally, the behavior of a function near a critical point can be characterized in terms of the eigenvalues and eigenvectors of the Hessian matrix.  Let {(Subscript[v, i ])} represent the set of orthonormal eigenvectors and {(Subscript[λ, i ])} the corresponding eigenvalues.  Furthermore, let (Subscript[u, i]) = (Subscript[v, i]) .  (x - (Subscript[x, 0])) represent the displacements from the critical point with respect to the local coordinate system specified by the eigendirections.  The variation of f[x]near the critical point then separates

f[x] ≈ 
      f[(Subscript[x, 0])] + 
        (1/2) (∑_(i = 1)^n) (Subscript[λ, i]) ((Subscript[u, i])^2)  +  
        ⋯

into a sum of independent contributions.  If all eigenvalues have a common sign, the critical point is a local minimum if that sign is positive or a local maximum if negative, but if the signs are mixed the critical point is a saddle with respect to pairs of directions whose eigenvalues have opposite signs.  A vanishing eigenvalue indicates the existence of a ridge along the corresponding eigendirection.

This technique finds widespread applications in many fields.  For example, data analysis using least-squares fitting of a model function to data is an optimization problem in which one searches a multidimensional parameter space for the minimum squared deviation between the model and the data.  Alternatively, a linear-stability analysis of the behavior of systems obeying nonlinear differential equations near a fixed point is performed in this manner.

Write a gradient function using Map and a Hessian function using Outer.  [Do not use Table!]

Evaluate the eigenvalues and eigenvectors of the Hessian matrices for the examples shown above.

Each of the following functions has a critical point at the origin.  Diagonalize the Hessian matrix and classify the critical point as a minimum, maximum, saddle, or ridge.  Also, determine the angle between the eigensystem and the original coordinate axes.
a) x y Cos[x y]	
b) (((x - 2  y))^2) Exp[-((x^2) +  (y^2))]
c) (2  (x^2)  + (y^2)) Exp[-(2  (x^2) +  (y^2))]	
d) Cos[x +  y] Exp[-((x^2) +  (y^2))]

Analyze the behavior of the function f=(x+2y-3)(y-z+1)(x+z) near the origin by locating a critical point and constructing the eigensystem for the corresponding Hessian matrix.  [Hint: if you use Solve for g==0, inverse functions are returned because there exists a ridge for which the gradient equations are degenerate.]

Taylor series

Power series in one variable

The Taylor series expansion of a function f(x) about the point x = (Subscript[x, 0]) is produced using Series[f,{x,(Subscript[x, 0]),n}] where n is the highest power retained in the expansion.  Series can expand about ∞, will usually detect essential singularities, and will give a finite number of negative powers or logarithmic terms when necessary.

For example, the Lenard-Jones potential

Symbolize[(Subscript[v, 0])]; Symbolize[(Subscript[r, 0])]

v = -(Subscript[v, 0]) (2 (((((Subscript[r, 0])/r)))^6) - (((((Subscript[r, 0])/r)))^12));

is often used as a simple model of the interaction between two atoms in a diatomic molecule.  This parametrization obviously has a minimum at (Subscript[r, 0]).  Thus, it is useful to construct a power-series expansion that approximates the potential near its minimum.

vSeries = Series[v, {r, (Subscript[r, 0]), 4}]

Although this result appears, at first glance, to be printed as an ordinary expression, the presence of the final term representing all omitted terms of order n + 1 and higher shows that this expression is not simply an algebraic expression.  (How would one evaluate the truncation term?).  In fact, Series actually returns a SeriesData object whose arguments are the expansion variable, the expansion point, a list of terms, and a description of the omitted terms.

FullForm[vSeries]

A SeriesData object is converted to an ordinary algebraic expression using Normal[SeriesDataObject].

Normal[vSeries]

FullForm[Normal[vSeries]]

The coefficient of the term with order n is selected by SeriesCoefficient[series,n].  Thus, in this case the coefficient of (((r - (Subscript[r, 0])))^3) is

SeriesCoefficient[vSeries, 3]

Let U[y]=Tanh[y] or U[τ]=Tanh[1/τ] be two representations of the same function.  Compare the Taylor expansion of U[y] for small y with the Taylor expansion of U[τ] for large τ (expand about ∞).

Use series expansions to characterize the behavior of f[z] = 1 - (z/(Sqrt[(z^2) + (R^2)])) for both small and large values of z, assuming that R > 0.  Compare these approximations with the original function graphically and describe their regions of validity.  How does one obtain a useful approximation for z ≈ R?

Singularities

Mathematica will often be able to provide useful expansions near singularities.  Consider the familiar function Γ[x] with singularities for negative integers.

Gamma[-5]

If a function f[x] has a simple pole at x = (Subscript[x, 0]), the power series takes the form

f[x] ≈  
      (R/(x - (Subscript[x, 0]))) +  
        (Subscript[a, 0]) +  (Subscript[a, 1]) (x - 
              (Subscript[x, 0]))  +  ⋯

where R is the residue of f at (Subscript[x, 0]).  Thus, for this function we find

Series[Gamma[x], {x, -5, 2}]

Residue[Gamma[x], {x, -5}]

Series sometimes returns an expansion in fractional powers

Series[Sin[(Sqrt[x])], {x, 0, 2}]

or will isolate other characteristic behaviors.

Series[PolyLog[2, z], {z, 1, 2}] // Normal // Simplify

Useful information can often be obtained even near an essential singularity.

Series[Exp[Sin[x]], {x, ∞, 2}]

Limit[Exp[Sin[x]], x->∞]

Here we are told that even though this function oscillates wildly near its essential singularity, the range of oscillation is confined to a finite interval.

Directional limits can be obtained for some functions.

Limit[(1/x), x -> 0, Direction -> -1]

Limit[(1/x), x -> 0, Direction -> +1]

Find the residue of x tan(x) at (π/2).

Evaluate the limit of Tanh[(1/x)] as x->0.  How does the result depend upon direction?

Compare the behaviors of BesselJ[2,z] and BesselK[2,z] as z->∞.

Multivariate power series

Series[f, {x, (Subscript["\!\(\*StyleBox[\"\\\"x\\\"\", \"TI\"]\)", 0]), (Subscript["\!\(\*StyleBox[\"\\\"n\\\"\", \"TI\"]\)", x])}, {y, (Subscript["\!\(\*StyleBox[\"\\\"y\\\"\", \"TI\"]\)", 0]), (Subscript["\!\(\*StyleBox[\"\\\"n\\\"\", \"TI\"]\)", y])} ] expands a function of two variables with respect to x first, then with respect to y.  Additional variables can be added to the sequence also.  Note that the reference guide describes the order of these operations incorrectly; the net result is usually the same in either order, but the presentation depends upon order.

Series[Log[1 - x y], {x, 1, 3}]

Series[%, {y, 0, 3}]

xySeries  = Series[Log[1 - x y], {x, 1, 3}, {y, 0, 3}]

yxSeries = Series[Log[1 - x y], {y, 0, 3}, {x, 1, 3}]

Normal[xySeries] - Normal[yxSeries] // Simplify

Symbolic integration

Symbolic integration in one dimension

The basic syntax for symbolic one-dimensional integration is Integrate[f[x],x] for an indefinite integral or Integrate[f[x],{x,a,b}] for a definite integral over the range {a, b}.  The traditional forms can be entered from the BasicInput palette as

∫f[x] dx

for indefinite, or

(∫_a^b) f[
          x] dx

for definite integrals, where d is not an ordinary d but a differential operator that can be entered as EscddEsc.  Most integrals one might find in standard tables will be evaluated successfully, although it is usually prudent to simplify integrands before requesting integration.

Indefinite integrals

Notice that indefinite integrals

∫(x^n) (a^x) \[DifferentialD]x

are returned without constants of integration and can be checked easily by differentiation

D[%, x] // Simplify

although simplification is often needed to recover the original form.  Also notice that Mathematica makes no a priori assumptions about the parameters involved in indefinite integrals

∫Exp[-a x] Cos[b x] \[DifferentialD]x  // 
      Simplify

so that the results apply to real or complex parameters.  On the other hand, it does not check special cases either.  Thus, the following integral

∫(x^μ) \[DifferentialD]x

does not appear to be correct if the parameter μ happens to be -1.  Mathematica assumes that parameters in indefinite integrals should be interpreted generically.

Consider the generic indefinite integral
 (∫_a^b) (x^μ) dx=(((b^(μ + 1)) - (a^(μ + 1)))/(μ + 1))  
 for b > a > 0.  Show that the correct result for μ = -1 can be obtained by series expansion with respect to μ.

Evaluate and check the following integrals.
a) ∫((x^m)/(a + 
                (x^n))) dx 		
b) ∫ BesselJ[n,x] BesselJ[m,x] x dx

Definite integrals

Definite integrals with symbolic parameters in the integrand often return a conditional statement rather than a unique result.  Thus,

(∫_0^1) (x^n) \[DifferentialD]x

correctly restricts its evaluation to exponents greater than -1.   Similarly,

(∫_0^∞) Exp[a x] \[DifferentialD]x

can be evaluated only when the parameter a has negative real part.  If we happen to know that the parameter will always satisfy this condition, we can incorporate this assumption in the integration command using the option Assumptions as follows.

Integrate[Exp[a x], {x, 0, ∞}, 
      Assumptions -> Re[a] < 0]

Assumptions can be supplied to typeset integrals using Assuming[assumptions,expression] where assumptions is a list of assumptions to be used in evaluation of an expression.  Thus, the preceding result can also be obtained as follows.

Assuming[
      Re[a] < 0, (∫_0^∞) Exp[
            a x] \[DifferentialD]x]

Alternatively, a reckless person could disable conditional-testing using

Integrate[Exp[a x], {x, 0, ∞}, 
      GenerateConditions -> False]

but that practice is apt to generate nonsense and is strongly discouraged unless you are certain that the conditions would be satisfied for your parameters.

Definite integrals with several parameters are often expressed in terms of composite conditions:

result = 
      Integrate[
        Exp[-a x]/((b^2) + (x^2)), {x, 0, ∞}]

Assuming that the conditions are satisfied, the desired integral can be extracted from the second part of result.

result[[2]]

In previous versions definite integrals with symbolic limits were treated in much the same manner as indefinite integrals and, hence, were evaluated without regard to restrictions upon the values of their parameters or limits.  In version 5 we obtain a compilation of possibilities according to the nature of the symbolic limits.

(∫_a^b) (x^n) \[DifferentialD]x

However, this procedure does not appear to examine the nature of the parameters in the integrand and is not correct for n = -1 even when the limits are real.

Compare
a) (∫_0^π)((Cos[t])^((3/2)))dt
b) (∫_0^π)((Cos[t])^(n + (1/2)))dt

Evaluate the following integrals.
a) (∫_0^∞)(( 
                  x BesselJ[n, 
                    x])/(((1 + (x^2)))^(m + (1/2)))) dx

Evaluate the following integrals.
a) (∫_0^(2 π))Cos[
              z Sin[t]] dt	
b) (∫_0^(2 π))Exp[
              z Sin[t]] dt

Evaluate \(TraditionalForm\`∫\_0\%ωD\( ω\^3\/\(e\^\(ω/T\) - 1\)\) dω\) which occurs in the Debye theory of lattice vibrations.

Evaluate \(∫\_0\%∞\(\( ω\^4\) e\^ω\)\/\((e\^ω - 1)\)\^2\).

rms radius for spherical distribution (symbolic)

\(TraditionalForm\`rms = \@\(\(∫\_0\%∞ ρ\ \(r\^4\) dr\)\/\(∫\_0\%∞ ρ\ \(r\^2\) dr\)\)\)

Principal value

In general one must treat singularities of an integrand with great care.  For the present purposes it suffices to mention that Mathematica provides an option for the Cauchy principal value in which a small interval surrounding a singularity is excised and the integral is evaluated in the limit that the excluded interval shrinks to zero.

pv = 
      Integrate[
        Exp[x]/((x + 1) (x + 2)), {x, -∞, 0}, 
        PrincipalValue -> True]

pv // N

Line integrals

A line integral evaluates and integrates a function along a specified path through a multidimensional space, returning a scalar result.  Here we consider two types of line integral.  The integral of scalar function against the path length for a curve in two dimensions takes the form

(∫_C) f[x, 
          y] ds

where C denotes a curve and \(ds = \@\(\((dx)\)\^2 + \((dy)\)\^2\)\) represents the path length along the curve.  The curve is generally represented by the parametric equations C = {x -> x[t], y -> y[t]} such that the path length becomes \(TraditionalForm\`ds = dt \@\(\((dx\/dt)\)\^2 + \((dy\/dt)\)\^2\)\).  For example, if f represents mass density, an integral of this type could be used to determine the total mass of a rope with nonuniform density or cross-sectional area.

Similarly, the component of a vector field that is parallel to a curve is integrated using

(∫_C) (Overscript[f, ⇀])[x, 
          y] . d(Overscript[s, ⇀])

where \(TraditionalForm\`ds\&⇀ = dt {dx\/dt, dy\/dt}\) is the differential displacement along the curve.  For example, the work performed by a force (Overscript[f, ⇀]) that acts on a particle whose trajectory is described by the parametric curve C = {x -> x[t], y -> y[t]} is given by an integral of this type.
The following function evaluates these types of path integrals symbolically.  Here we present a relatively simple implementation and do not attempt to cover all possible methods for specifying the function or path.  The path-integration function can be generalized later if needed.  We assume that the function is an implicit function of the two variables represented by patterns x_ and y_.  The path is defined using a pattern argument path:{pattern rules} in which the pattern rules are a list of rewrite rules of the form x_->u_ that replace the variable on the lhs by the pattern on the rhs; hence, blank patterns are used on both sides of the arrow.  The parameter defining the progress along the curve is given by t_.  Two versions are given, one for scalar functions and the other for vector functions in two dimensions.  Obviously, it is simple matter to generalize to arbitrary dimensionality.

ClearAll[pathIntegrate]

\(\(pathIntegrate::usage = "pathIntegrate[f,{x->u,y->v},t,tmin,tmax,opts] evaluates the line integral of the scalar function f along a curve expressed as a parametric function of t. 

 pathIntegrate[{fx,fy},{x->u,y->v},t,tmin,tmax] evaluates the line integral of the tangential component of a vector field along a curve expressed as a parametric function of t.";\)
		
		\)
\(\(\(pathIntegrate[f_ /; \((Not[VectorQ[f]])\), 
          path : {x_ -> u_, y_ -> v_}, t_, tmin_, tmax_, 
          opts___Rule]\  := 
        Integrate[
          f\ √\((D[u, t]\^2 + D[v, t]\^2)\) /. path, {t, tmin, tmax}, 
          opts];\)\(
\)
    \)\)
pathIntegrate[field : {f_, g_}, path : {x_ -> u_, y_ -> v_}, 
        t_, tmin_, tmax_, opts___Rule]  := 
      Integrate[(f D[u, t] + g D[v, t]) /. path, {t, tmin, tmax}, 
        opts]

Examples

pathIntegrate[(1/((x^2) + (y^2))), {x -> Cos[t], y -> Sin[t]}, 
      t, 0, 2  π]

pathIntegrate[(1/((x^2) + (y^2))), {x -> t, y -> t}, t, 1, 
      2]

pathIntegrate[(1/((x^2) + (y^2))), {x -> a  + t, y -> b t}, 
      t, 0, 1, Assumptions -> {a ∈ Reals, 
          b ∈ Reals}]

pathIntegrate[{(x^2), y}, {x -> t, y -> t}, t, 0, 
      2]

pathIntegrate[{(x^2), y}, {x -> 2  Sin[t], y -> 2  Sin[t]}, t, 
      0, (π/2)]

pathIntegrate[{y, x}, {x -> a (t^n), y -> b (t^m)}, t, 0, 
      1]

Write line integration functions applicable to arbitrary dimensionality.

Proper time

Multiple integrals

Multiple integrals are entered simply by replacing the integration variable (and limits) with a sequence of integration variables.  The integrals are evaluated from right to left, so that the integration limits for earlier variables can be functions of the later variables to define a region of integration that is not simply rectangular.  Thus,

Clear[f, x, y, xmin, xmax, ymin, ymax];
Integrate[f[x,  y], {y, ymin, ymax}, {x, xmin, xmax}]

A trivial indefinite double integral is

Integrate[x (y^2), x, y]

while the integration over a rectangular region is given by

Integrate[x (y^2), {x, xmin, xmax}, {y, ymin, ymax}]

The following example integrates within the triangular wedge defined by 0 ≤ x ≤ y and 0 ≤ y ≤ 2.  Thus,

Integrate[x (y^2), {y, 0, 2}, {x, 0, y}]

(∫_0^2)(∫_0^y) 
          x  (y^2) \[DifferentialD]x \[DifferentialD]y

A more interesting example is the function ylmCoeff which uses the orthonormality of spherical harmonics to produce expansion coefficients for angular functions.

ylmCoeff[λ_, m_, f_] :=  (((-1))^m) 
        Integrate[
          f SphericalHarmonicY[λ, -m, θ, ϕ] 
            Sin[θ], {θ, 0, π}, {ϕ, 0, 
            2  π}]

Table[
      ylmCoeff[2, m, 
        SphericalHarmonicY[2, -1, θ, ϕ]], {m, -2, 
        2}]

Table[
        ylmCoeff[2, m, Cos[θ]^2 Sin[2  ϕ]], {m, -2, 2}] // 
      Simplify

Table[ylmCoeff[2, m, Cos[θ] Sin[ϕ]], {m, -2, 2}] // 
      Simplify

Area using change of variables

Contour integration

Contour integrals are very similar to line integrals except that the path lives in the complex plane.  If z is a complex variable, a contour can be defined by a parametric function z -> ζ[t] where t varies between values that take ζ[t] from the beginning to the end of the desired path.  The following function replaces the complex variable by its value on the contour and converts the displacement from dz to dt by means of the partial derivative D[ζ, t].

contourIntegrate[f_, 
        contour : (z_ -> ζ_), {t_, tmin_, tmax_}] := 
      Integrate[f D[ζ, t] /. contour, {t, tmin, tmax}]

The theorem of residues tells us that the integral of a function that is regular (analytic and single-valued) on a simple, closed, counterclockwise contour and the enclosed domain, except for a finite number of isolated singularities, is equal to 2  π itimes the sum of the residues for the singularities within the contour.  Hence, we test our function using a function whose one simple pole has unit residue

contourIntegrate[(1/(z - (1/2))), 
      z -> Cos[t] + I Sin[t], {t, 0, 2  π}]

and obtain the correct answer.  Integration along the more complicated limaçon contour suggested by Gass should give 4  π  i because the contour encircles the pole twice.

plotContour[f_, {t_, tmin_, tmax_}, opts___Rule] := 
      ParametricPlot[Evaluate[{Re[f], Im[f]}], {t, tmin, tmax}, 
        opts]

Show[{plotContour[Exp[I t] (2  Cos[t] + 1), {t, 0, 2  π}, 
            AspectRatio -> Automatic, DisplayFunction -> Identity], 
          Graphics[{AbsolutePointSize[10], Point[{(1/2), 0}]}]}, 
        DisplayFunction -> $DisplayFunction];

contourIntegrate[(1/(z - (1/2))), 
        z -> Exp[I t] (2  Cos[t] + 1), {t, 0, 
          2  π}] // Simplify

Version 3 gave only half the correct answer, Version 4 gave the correct answer, but Version 5.1 again gives an incorrect answer, this time including a spurious real part.  One must always adhere to the nuclear arms control adage:

Evaluate a contour integral in which a simple pole is enclosed by a circle that is traversed twice.  Compare with the situation in which two loops of different sizes enclose the same pole.

Integrate (1/(z - (1/2))) - (2/(z - i - 1)) on both of the following contours: z->Cos[t]+i Sin[t] and z->Exp[i t](2Cos[t]+1) for {t,0,2π}.  Comment on your results.

Numerical integration

Basic syntax

Although Mathematica will often succeed in producing symbolic integrals for some quite unpromising integrands, many interesting integrals cannot be performed symbolically.  Furthermore, it is often more efficient to evaluate an integral numerically even when a symbolic result is available because many symbolic integrals are expressed in terms of rather complicated and unfamiliar functions anyway.  Not surprisingly, the function Mathematica provides for numerical integration is called NIntegrate and is used with syntax virtually identical to that for symbolic integration.  Of course, the integrand must evaluate to a number, which sometimes requires wrapping the function Evaluate around the integrand.

Integrate[Exp[ Cos[t]^2], {t, 0, 2  π}]

% // N

NIntegrate[Exp[ Cos[t]^2], {t, 0, 2  π}]

Here is an example that Version 4 can integrate symbolically, but Version 3 cannot.

integral = 
      Integrate[(((x^2) Exp[-x])/(1 + (x^3))), {x, 
          0, ∞}]

However, few would find the symbolic result helpful!  Furthermore, numerical evaluation of this symbolic results does take measurable time

{time1, ans1} = Timing[N[integral]]

while direct numerical evaluation of the integral is faster

{time2, ans2} = 
      Timing[NIntegrate[(((x^2) Exp[-x])/(1 + (x^3))), {x, 
            0, ∞}]]

and gives the same result to high accuracy

ans1 - ans2

Although this difference in speed is unimportant for a single evaluation, one often needs to evaluate an integral many times within a larger program or loop for which efficiency does become important.

Finally, one should always remember that no numerical algorithm is guaranteed to succeed for any input.  The integrand for the present example is so well behaved that numerical integration is obviously straightforward

Plot[(((x^2) Exp[-x])/(1 + (x^3))), {x, 0, 5}];

but Murphy's Law decrees that the integral you need for your work will not be as compliant.  You must treat any singularities or unusually rapid changes carefully.  Numerical integration often has difficulties with infinite ranges if the integrand does not decrease sufficiently rapidly or oscillates.  Often one can improve the convergence by a suitable change of variable or by subdividing the integration range.  NIntegrate also provides several options

Options[NIntegrate]

which can often tame a slowly convergent integral by brute force, but insight into the properties of your integral is usually preferred.  However, a systematic treatment of reluctant integrals is beyond the scope of this course.

Compare symbolic and numerical evaluations of \(∫\_0\%1\(\( ω\^4\) e\^ω\)\/\((e\^ω - 1)\)\^2\).

Evaluate \(∫\_0\%∞\(\( ω\^4\) e\^ω\)\/\((e\^ω - 1)\)\^2\) numerically.  [Hint: use LogPlot to examine the integrand.]

Example: finite Fourier-Bessel transform

Consider a spherically symmetric three-dimensional density of charge or mass that is relatively large near the origin and vanishingly small for large distances.  An example dear to the author's heart is the density of atomic nuclei, which is reasonably well approximated by the simple function

Clear[ρ];
ρ =  ((1 + w (((r/c))^2))/(1 + Exp[(r - c)/a]));

A typical case is shown below.  Note that the central density is normalized artificially to unity.  The density is fairly constant for small radii and then decays to zero in a distance of order a near the surface radius c.

values  =  {c -> 2.5, a -> 0.5, 
          w -> 0.2};

Plot[Evaluate[ρ /. values], {r, 0, 8}];

When analyzing measurements of such a density distribution, one commonly employs a Fourier-Bessel transform

\(TraditionalForm\`fbe[
        r]\  = \ \(\(∑\+\(n = 1\)\%nmax\( 
                a\_n\) 𝒿\_0[\(n\ π\ r\)\/R]\ \ for\ r\)\(≤\)\(R\)\(\ \)\)\)

within a limited volume whose radius is several times the root-mean-square radius of the distribution.  Using known properties of the spherical Bessel functions, one can show that the expansion coefficients are given by the integral

\(TraditionalForm\`a\_n = \(2\/\(\(R\^3\) 𝒿\_1[\(n\ π\ r\)\/R]\^2\)\) \(∫\_0\%R ρ[
              r] 𝒿\_0[\(n\ π\ r\)\/R] \(r\^2\) dr\)\)

Thus, we wish to construct a function which evaluates the expansion coefficients by numerical integration.  We also test this function by comparing the expansion with the original function.

Clear[sbessj, fbeCoeff];
sbessj[n_, z_  /; z ≠ 0]  := 
        Sqrt[ Pi/(2  z)] BesselJ[n + (1/2), z];
sbessj[0, 0]  :=  1;
sbessj[n_ /; n > 0, 0 ]  :=  0

fbeCoeff[ρ_, params_, n_, 
        R_] :=   (2/((R^3) sbessj[1, n π ]^2)) 
        NIntegrate[
          Evaluate[
            (r^2) ρ sbessj[0, ((n π r)/R)] /. params], {r, 0, 
            R}]

fbe[r_, coeffs_, R_] :=  
      coeffs . Table[
          sbessj[0, ((n π r)/R)], {n, 1, Length[coeffs]}]

coeffs = 
      Table[fbeCoeff[ρ, {c -> 2.5, a -> 0.5, w -> 0.2}, 
          n, 8. ], {n, 1, 10}]

Plot[
        Evaluate[(ρ - fbe[r, coeffs, 8. ]) /. {c -> 2.5, 
              a -> 0.5, w -> 0.2}], {r, 0, 8}];

For this choice of cut-off radius we achieve better than 1% accuracy with 10 terms of the Fourier-Bessel expansion.

rms radius for spherical distribution (numerical)

\(TraditionalForm\`rms = \@\(\(∫\_0\%∞ ρ\ \(r\^4\) dr\)\/\(∫\_0\%∞ ρ\ \(r\^2\) dr\)\)\)

Numerical contour integration

Contour integration along a path defined as a series of connected line segments can be performed using the syntax NIntegrate[f, {z, (Subscript[z, 1]), (Subscript[z, 2]), ⋯, (Subscript[z, n])}] where the (Subscript[z, i]) are complex numbers that plot the contour.  For example, the following integral is performed along a square contour centered on the origin with edges parallel to the real and imaginary axis.

NIntegrate[
      Exp[z]/z, {z, -1 - I, 1 - I, 
        1 + I, -1 + I, -1 - I}]

Because we chose a function exhibiting only one simple pole with unit residue inside the contour, the result is just 2  π i.  Complex analysis tells us that the same result should be obtained for any other contour which encompasses the same poles and which can be deformed continuously into the first without crossing singularities or branch cuts.  Thus, the diamond-shaped contour

NIntegrate[
      Exp[z]/z, {z, -1, -I, 
        1, I, -1}]

does indeed give the same result.

More general contours can be defined in terms of the same type of replacement rules we employed above for symbolic contour integration.

contourNIntegrate[f_, 
        contour : (z_ -> ζ_), {t_, tmin_, tmax_}] := 
      NIntegrate[
        Evaluate[f D[ζ, t] /. contour], {t, tmin, tmax}]

Applying this over a circle,

contourNIntegrate[Exp[z]/z, 
        z -> Cos[t] + I Sin[t], {t, 0, 2  π}] // 
      Chop

we again obtain the same result when the contour encloses a single pole.  Returning to our earlier example in which the pole is circled twice, we now obtain the correct answer because numerical integration counts both loops.

contourNIntegrate[(1/(z - (1/2))), 
        z -> Exp[I t] (2  Cos[t] + 1), {t, 0, 
          2  π}] // Chop

Nor is the spurious real part produced by the symbolic integral present in the numerical integral.

Numerically integrate (1/(z - (1/2))) - (2/(z - i - 1)) on both of the following contours for {t,0,2π}: 
a) z->Cos[t]+i Sin[t]  		
[Hint: do not separate the integrand.]
b) z->Exp[i t](2Cos[t]+1)	
[Hint: you may need to handle the two terms separately for the two-loop contour because nearby poles tend to cause problems with numerical stability.]  
Comment on your results.

\(TraditionalForm\`f[n, α] = 
      Sin[n\ α]\^2\/Sin[α]\^2\)

α  =  (d/λ) 
        Sin[θ]

\(TraditionalForm\`m\ ∂\^2 r\/∂t\^2 == \(-\(∂V\/∂r\)\)\)

V[r] ≈ 
      (Subscript[V, 0]) + (1/2) 
          k (((r - (Subscript[r, 0])))^2) + ⋯

V[
        r] = -(( d)/(1 - 
                2  (c^2))) (Exp[-2 a(r - b)/c]  - 
            2  (c^2) Exp[-a(r - b)/c])

\(TraditionalForm\`ϕ[r\&⇀] - ϕ[
          r\&⇀\_0] = \(-\(∫\_\(r\&⇀\_0\)\%\(r\&⇀\)f\&⇀ . dr\&⇀\)\)\)