Chapter 3

Calculus

Abstract

Chapter 3 introduces Mathematica's calculus capabilities. The examples used to illustrate the various functions are similar to examples typically seen in a traditional calculus sequence. If you have trouble typing commands correctly, use the buttons on the Basic Math Input or Basic Math Assistant palettes to help you create templates in standard mathematical notation that you can evaluate.

Keywords

Limits; Differential calculus; Integral calculus; Series calculus; Multi-variable calculus

3.1 Limits and Continuity

One of the first topics discussed in calculus is that of limits. Mathematica can be used to investigate limits graphically and numerically. In addition, the Mathematica command Limit[f[x],x->a] attempts to compute the limit of y=f(x) as x approaches a, limx→a⁡f(x), where a can be a finite number, ∞ (Infinity), or −∞ (-Infinity). The arrow “->” is obtained by typing a minus sign “-” followed by a greater than sign “>”.

Remark 3.1

To define a function of a single variable, f(x)=expressioninx, enter f[x_]=expression in x. To generate a basic plot of y=f(x) for a⩽x⩽b, enter Plot[f[x],{x,a,b}].

3.1.1 Using Graphs and Tables to Predict Limits

Example 3.1

Use a graph and table of values to investigate limx→0⁡sin⁡3xx.

Solution

We clear all prior definitions of f, define f(x)=(sin⁡3x)/x, and then graph y=f(x) on the interval [−π,π] with Plot.

Clear[f] clears all prior definitions of f, if any. Clearing function definitions before defining new ones helps eliminate any possible confusion and/or ambiguities.

Clear [ f ]

f [ x _ ] = Sin [ 3 x ] x ;

Plot [ f [ x ] , { x , − 2 π , 2 π } ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ { . 05 , . 26 , 1 , . 27 } ] } } ,

PlotRange → All ]

From the graph shown in Fig. 3.1, we might, correctly, conclude that limx→0⁡sin⁡3xx=3. Further evidence that limx→0⁡sin⁡3xx=3 can be obtained by computing the values of f(x) for values of x “near” x=0. In the following, we use RandomReal. RandomReal[{a,b}] returns a “random” real number between a and b. Because we are generating “random” numbers, your results will differ from those obtained here, to define xvals to be a table of 6 “random” real numbers. The first number in xvals is between −1 and 1, the second between −1/10 and 1/10, and so on.

Figure 3.1 Graph of f(x)=(sin⁡3x)/x on the interval [−2π,2π] in Wisconsin gold.

Remark 3.2

Throughout Mathematica by Example we illustrate how different options such as those that affect the coloring and labeling of a graphic are used, which include options such as PlotStyle, PlotLabel, and AxesLabel.

xvals=Table[RandomReal[{−10−n,10−n}],{n,0,5}]

{ − 0.710444 , − 0.0134926 , − 0.00161622 , − 0.000546629 , 0.0000262461 ,

3.6298466289834744 ` * − ∧ 6 }

We then use Map to compute the value of f(x) for each x in xvals. Map[f,{x1,x2,...,xn}] returns the list {f(x1),f(x2),...,f(xn)}. We use Table to display the results in tabular form. Generally, list[[i]] returns the ith element of list while

Table[f[i],{i,start,finish,stepsize}]

computes each value of f(i) from start to finish in increments of stepsize. To create a table consisting of n equally spaced values, use stepsize=(finish-start)/(n-1). TableForm attempts to display a table form in a standard format such as the row-and-column format that follows.

pairs = Table [ { xvals [ [ i ] ] , fvals [ [ i ] ] } , { i , 1 , 6 } ] ;

TableForm [ pairs ]

− 0.710444 1.19217 − 0.0134926 2.99918 − 0.00161622 2.99999 − 0.000546629 3 . 0.0000262461 3 . 3.6298466289834744 ` * − ∧ 6 3 .

From these values, we might again correctly deduce that limx→0⁡sin⁡3xx=3. Of course, these results do not prove that limx→0⁡sin⁡3xx=3 but they are helpful in convincing us that limx→0⁡sin⁡3xx=3. □

For piecewise-defined functions, you can either use Mathematica's conditional command (/;) to define the piecewise-defined function or use Piecewise.

Example 3.2

If h(x)={x2+x, if x⩽01−x, if 0<x<32x2−15x+25, if 3⩽x⩽515−2x, if x>5, compute the following limits: (a) limx→0⁡h(x), (b) limx→3⁡h(x), (c) limx→5⁡h(x).

Solution

We use Mathematica's conditional command, /; to define h. We must use delayed evaluation (:=) because h(x) cannot be computed unless Mathematica is given a particular value of x. The first line of the following defines h(x) to be x2+x for x⩽0, the second line defines h(x) to be 1−x for 0<x<3 and so on. Notice that Mathematica accidentally connects (0,0) to (0,1) and then (5,0) to (5,5). (See Fig. 3.2 (a).) The delayed evaluation is also incompatible with Mathematica's Limit function.

Figure 3.2 (a) Plot does not catch the breaks in the piecewise defined function (Badger red). (b) If you use Piecewise, Plot can catch jumps (Illinois blue).

The plots p1 and p2 are not displayed because a semi-colon is included at the end of each Plot command.

Clear [ h ]

h [ x _ ] := x ∧ 2 + x /; x ⩽ 0

h [ x _ ] := 1 − x /; 0 < x < 3

h [ x _ ] := 2 x ∧ 2 − 15 x + 25 /; 3 ⩽ x ⩽ 5

h [ x _ ] := 15 − 2 x /; x > 5

p1 = Plot [ h [ x ] , { x , 0 , 6 } ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ { . 002 , 1 , . 85 , . 06 } ] } } ]

To avoid these problems and see the discontinuities in our plots, we redefine h using Mathematica's Piecewise function as follows.

Clear [ h ]

h [ x _ ] := Piecewise [ { { x ∧ 2 + x , x < 0 } , { 1 − x , 0 < x < 3 } , { 2 x ∧ 2 − 15 x + 25 , 3 < = x ⩽ 5 } ,

{ 15 − 2 x , x > 5 } } ] ;

p2 = Plot [ h [ x ] , { x , − 2 , 6 } , PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 1 , . 86 , . 24 , . 09 ] } } ,

PlotRange→All];

Show [ GraphicsRow [ { p1 , p2 } ] ]

Notice that when we execute the Plot command, Mathematica “catches” the breaks between (0,0) and (0,1) and then at (5,0) and (5,5) shown in Fig. 3.2 (b).

From Fig. 3.2, we see that limx→0⁡h(x) does not exist, limx→3⁡h(x)=−2, and limx→5⁡h(x) does not exist. □

When limits exist, you can often use Limit[f[x],x->a] (where a may be ±Infinity) to compute limx→a⁡f(x). Thus, for the previous example we see that

Limit [ h [ x ] , x → 3 ]

−2

is correct. On the other hand

Limit [ h [ x ] , x → 5 ]

5

is incorrect. We check by computing the right hand limit, limx→5+⁡h(x), using the Direction->-1 option in the Limit command and then the left limit, limx→5−⁡h(x), using the Direction->1 in the Limit command.

Limit [ h [ x ] , x → 5 , Direction → 1 ]

0

Limit [ h [ x ] , x → 5 , Direction → − 1 ]

5

We follow the same procedure for x=0

Limit [ h [ x ] , x → 0 ]

1

Limit [ h [ x ] , x → 0 , Direction → 1 ]

0

Limit [ h [ x ] , x → 0 , Direction → − 1 ]

1

3.1.2 Computing Limits

Some limits involving rational functions can be computed by factoring the numerator and denominator.

Example 3.3

Compute limx→−9/2⁡2x2+25x+7272−47x−14x2.

Solution

We define frac1 to be the rational expression 2x2+25x+7272−47x−14x2. We then attempt to compute the value of frac1 if x=−9/2 by using ReplaceAll (/.) to evaluate frac1 if x=−9/2 but see that it is undefined.

Factoring the numerator and denominator with Factor, Numerator, and Denominator, we see that

lim x → − 9 / 2 ⁡ 2 x 2 + 25 x + 72 72 − 47 x − 14 x 2 = lim x → − 9 / 2 ⁡ ( x + 8 ) ( 2 x + 9 ) ( 8 − 7 x ) ( 2 x + 9 ) = lim x → − 9 / 2 ⁡ x + 8 8 − 7 x .

The fraction (x+8)/(8−7x) is named frac2 and the limit is evaluated by computing the value of frac2 if x=−9/2

Factor [ Numerator [ frac1 ] ]

Factor [ Denominator [ frac1 ] ]

( 8 + x ) ( 9 + 2 x )

− ( 9 + 2 x ) ( − 8 + 7 x )

frac2 = Cancel [ frac1 ]

− 8 − x − 8 + 7 x

frac2 /. x → − 9 2

7 79

or by using the Limit function on the original fraction.

Limit [ frac1 , x → − 9 / 2 ]

7 79

We conclude that

lim x → − 9 / 2 ⁡ 2 x 2 + 25 x + 72 72 − 47 x − 14 x 2 = 7 79 .

 □

As stated previously, Limit[f[x],x->a] attempts to compute limx→a⁡f(x),

Limit[f[x],x->a,Direction->1] attempts to compute limx→a−⁡f(x), and

Limit[f[x],x->a,Direction->-1] attempts to compute limx→a+⁡f(x). Generally, a can be a number, ±Infinity (±∞), or another symbol.

Thus, entering

Limit [ 2 x 2 + 25 x + 72 72 − 47 x − 14 x 2 , x → − 9 2 ]

7 79

computes limx→−9/2⁡2x2+25x+7272−47x−14x2=7/79.

Example 3.4

Calculate each limit: (a) limx→−5/3⁡3x2−7x−2021x2+14x−35; (b) limx→0⁡sin⁡xx; (c) limx→∞⁡(1+zx)x; (d) limx→0⁡e3x−1x; (e) limx→∞⁡e−2xx; and (f) limx→1+⁡(1ln⁡x−1x−1).

Solution

In each case, we use Limit to evaluate the indicated limit. Entering

Limit [ 3 x 2 − 7 x − 20 21 x 2 + 14 x − 35 , x → − 5 3 ]

17 56

computes limx→−5/3⁡3x2−7x−2021x2+14x−35=1756; and entering

Limit [ Sin [ x ] x , x → 0 ]

1

computes limx→0⁡sin⁡xx=1. Mathematica represents ∞ by Infinity. Thus, entering

Limit [ ( 1 + z / x ) ∧ x , x → Infinity ]

e z

computes limx→∞⁡(1+zx)x=ez. Entering

Limit [ ( Exp [ 3 x ] − 1 ) / x , x → 0 ]

3

computes limx→0⁡e3x−1x=3. Entering

Limit [ Exp [ − 2 x ] Sqrt [ x ] , x → Infinity ]

0

computes limx→∞⁡e−2xx=0, and entering

Because ln⁡x is undefined for x⩽0, a right-hand limit is mathematically necessary, even though Mathematica's Limit function computes the limit correctly without the distinction.

Limit [ 1 / Log [ x ] − 1 / ( x − 1 ) , x → 1 , Direction → − 1 ]

1 2

computes limx→1+⁡(1ln⁡x−1x−1)=12. □

3.1.3 One-Sided Limits

As illustrated previously, Mathematica can compute certain one-sided limits. The command Limit[f[x],x->a,Direction->1] attempts to compute limx→a−⁡f(x) while Limit[f[x],x->a,Direction->-1] attempts to compute limx→a+⁡f(x).

Example 3.5

Compute (a) limx→0+⁡|x|/x; (b) limx→0−⁡|x|/x; (c) limx→0+⁡e−1/x; and (d) limx→0−⁡e−1/x.

Solution

Even though limx→0⁡|x|/x does not exist, limx→0+⁡|x|/x=1 and limx→0−⁡|x|/x=−1, as we see using Limit together with the Direction->-1 and Direction->1 options, respectively.

Limit [ Abs [ x ] x , x → 0 , Direction → 1 ]

Limit [ Abs [ x ] x , x → 0 , Direction → − 1 ]

−1

1

The Direction->-1 and Direction->1 options are used to calculate the correct values for (c) and (d), respectively. For (c), we have:

Limit [ 1 x , x → 0 , Direction → − 1 ]

∞

Technically, limx→0⁡e−1/x does not exist (see Fig. 3.3 (a)) so the following is incorrect.

Figure 3.3 (a) Graph of y = e^−1/x (Army green). (b) Graph of y=e−1/x2 (Navy blue).

Limit [ Exp [ − 1 / x ] , x → 0 ]

0

However, using Limit together with the Direction option gives the correct left and right limits.

Limit [ Exp [ − 1 / x ] , x → 0 , Direction → 1 ]

∞

Limit [ Exp [ − 1 / x ] , x → 0 , Direction → − 1 ]

0

We confirm these results by graphing y=e−1/x with Plot in Fig. 3.3 (a). In (b), we also show the graph of y=e−1/x2 in Fig. 3.3 (b), which is further discussed in the exercises.

p1 = Plot [ Exp [ − 1 / x ] , { x , − 5 , 5 } , PlotStyle → CMYKColor [ . 64 , . 47 , 1 , . 4 ] ,

PlotLabel → “(a)” ] ;

p2 = Plot [ Exp [ − 1 / x ∧ 2 ] , { x , − 5 , 5 } ,

PlotStyle → CMYKColor [ 1 , . 38 , 0 , . 64 ] ,

PlotLabel → “(b)” ] ;

Show[GraphicsRow[{p1,p2}]] □

The Limit command together and its options (Direction->1 and Direction->-1) are “fragile” and should be used with caution because the results can be unpredictable. It is wise to check or confirm results using a different technique for nearly all problems encountered.

3.1.4 Continuity

Definition 3.1

The function y=f(x) is continuous at x=a if

1.  limx→a⁡f(x) exists;

2.  f(a) exists; and

3.  limx→a⁡f(x)=f(a).

Note that the third item in the definition means that both (1) and (2) are satisfied. But, if either of (1) or (2) is not satisfied the function is not continuous at the number in question. The function y=f(x) is continuous on the open interval I if f(x) is continuous at each number a contained in the interval I. Loosely speaking, the “standard” set of functions (polynomials, rational, trigonometric, etc...) are continuous on their domains.

Remark 3.3

Be careful with regard to this. For example, since limx→0−⁡x does not exist, many would say that f(x)=x is right continuous at x=0.

Example 3.6

For what value(s) of x, if any, are each of the following functions continuous? (a) f(x)=x3−8x; (b) f(x)=sin⁡2x; (c) f(x)=(x−1)/(x+1); (d) f(x)=(x−1)/(x+1).

Solution

(a) Polynomial functions are continuous for all real numbers. In interval notation, f(x) is continuous on (−∞,∞). (b) Because the sine function is continuous for all real numbers, f(x)=sin⁡2x is continuous for all real numbers. In interval notation, f(x) is continuous on (−∞,∞). (c) The rational function f(x)=(x−1)/(x+1) is continuous for all x≠−1. In interval notation, f(x) is continuous on (−∞,−1)∪(−1,∞). (d) f(x)=(x−1)/(x+1) is continuous if the radicand is nonnegative. In interval notation, f(x) is strictly continuous on (−∞,−1)∪(1,∞) but some might say that f(x) is continuous on (−∞,−1)∪[1,∞), where it is understood that f(x) is right continuous at x=1. We see this by graphing each function with the following commands. See Fig. 3.4. Note that in p3, the vertical line is not a part of the graph of the function—it is a vertical asymptote. If you were to redraw the figure by hand, the vertical line would not be a part of the graph.

Figure 3.4 (a) Polynomials (Michigan State green), (b) trigonometric (Ohio State gray), (c) rational (Alabama crimson), and (d) root functions (Auburn orange) are usually continuous on their domains.

p1 = Plot [ x ∧ 3 − 8 x , { x , − 5 , 5 } , PlotLabel → “(a)” ,

PlotStyle → CMYKColor [ . 82 , 0 , . 64 , . 7 ] ] ;

p2 = Plot [ Sin [ 2 x ] , { x , − 5 , 5 } , PlotLabel → “(b)” ,

PlotStyle → CMYKColor [ . 56 , . 47 , . 47 , . 15 ] ] ;

p3 = Plot [ ( ( x − 1 ) / ( x + 1 ) ) , { x , − 5 , 5 } , PlotLabel → “(c)” ,

PlotStyle → CMYKColor [ . 25 , . 97 , . 76 , . 18 ] ] ;

p4 = Plot [ Sqrt [ ( x − 1 ) / ( x + 1 ) ] , { x , − 5 , 5 } , PlotLabel → “(d)” ,

PlotStyle → CMYKColor [ 0 , . 61 , . 97 , 0 ] ] ;

Show[GraphicsGrid[{{p1,p2},{p3,p4}}]]  □

Computers are finite state machines so handling “interesting” functions can be problematic, especially when one must distinguish between rational and irrational numbers. We assume that if x=p/q is a rational number (p and q integers), p/q is a reduced fraction. One way of tackling these sorts of problems is to view rational numbers as ordered pairs, {a,b}. If a and b are integers, Mathematica automatically reduces a/b so Denominator[a/b] or a/b//Denominator returns the denominator of the reduced fraction; Numerator[a/b] or a/b//Numerator returns the numerator of the reduced fraction. If you want to see the points (x,f(x)) for which x is rational, we use ListPlot.

Example 3.7

Let f(x)={1/q, if x=p/q is rational0, if x is irrational. Create a representative graph of f(x).

Solution

You cannot see points: the measure of the rational numbers is 0, and the measure of the irrational numbers is the continuum, C. A true graph of f(x) would look like the graph of y=0. In the context of the example, we want to see how the graph of f(x) looks for rational values of x. We use a few points to illustrate the technique by using Table and Flatten to generate a set of ordered pairs.

Flatten[list,n] flattens list to level n.

In Mathematica, an ordered pair (a,b) is represented by {a,b}.

t1 = Flatten [ Table [ { n , m } , { n , 1 , 5 } , { m , 1 , 5 } ] , 1 ]

{ { 1 , 1 } , { 1 , 2 } , { 1 , 3 } , { 1 , 4 } , { 1 , 5 } , { 2 , 1 } , { 2 , 2 } ,

{ 2 , 3 } , { 2 , 4 } , { 2 , 5 } , { 3 , 1 } , { 3 , 2 } , { 3 , 3 } , { 3 , 4 } , { 3 , 5 } , { 4 , 1 } ,

{ 4 , 2 } , { 4 , 3 } , { 4 , 4 } , { 4 , 5 } , { 5 , 1 } ,

{ 5 , 2 } , { 5 , 3 } , { 5 , 4 } , { 5 , 5 } }

Next, we defined a function f. Assuming that a and b are integers, given an ordered pair {a,b}, f({a,b}) returns the point {a/b,1/(Reduced denominator of a/b)}

f [ { a _ , b _ } ] := { a / b , 1 / ( a / b // Denominator ) }

We use Map to compute the value of f for each ordered pair in t1. The resulting list is named t2.

t2 = Map [ f , t1 ]

{ { 1 , 1 } , { 1 2 , 1 2 } , { 1 3 , 1 3 } , { 1 4 , 1 4 } , { 1 5 , 1 5 } , { 2 , 1 } , { 1 , 1 } ,

{ 2 3 , 1 3 } , { 1 2 , 1 2 } , { 2 5 , 1 5 } , { 3 , 1 } , { 3 2 , 1 2 } , { 1 , 1 } ,

{ 3 4 , 1 4 } , { 3 5 , 1 5 } , { 4 , 1 } , { 2 , 1 } , { 4 3 , 1 3 } , { 1 , 1 } ,

{ 4 5 , 1 5 } , { 5 , 1 } , { 5 2 , 1 2 } , { 5 3 , 1 3 } , { 5 4 , 1 4 } , { 1 , 1 } }

Notice that t2 contains duplicate entries. We can remove them using Flatten but doing so does not affect the plot shown in Fig. 3.5 (a).

Figure 3.5 (a) After step 1 (Notre Dame gold). (b) After step 2 (USC cardinal). (c) Examining the numerator rather than the denominator (Oklahoma crimson). (d) The sine of the numerator (Texas burnt orange).

p1 = ListPlot [ t2 , PlotRange → { { 0 , 3 } , { 0 , 1 } } , AspectRatio → 1 ,

PlotStyle - > CMYKColor [ . 06 , . 27 , 1 , . 12 ] ,

PlotLabel → “(a)” ] ;

To generate a “prettier” plot, we repeat the procedure using more points. After entering each command the results are not displayed because we include a semicolon (;) at the end of each. See Fig. 3.5 (b).

t3 = Flatten [ Table [ { n , m } , { n , 1 , 300 } , { m , 1 , 200 } ] , 1 ] ;

t4 = Map [ f , t3 ] ;

p2 = ListPlot [ t4 , PlotRange → { { 0 , 3 } , { 0 , 1 } } , AspectRatio → 1 ,

PlotStyle - > CMYKColor [ . 07 , 1 , . 65 , . 32 ] ,

PlotLabel→“(b)”];

This function is interesting because it is continuous at the irrationals and discontinuous at the rationals.

We can consider other functions in similar contexts. In the following the y-coordinate is the numerator rather than the denominator. See Fig. 3.5 (c).

Clear [ f ]

f [ { a _ , b _ } ] := { a / b , a / b // Numerator } ;

t3 = Flatten [ Table [ { n , m } , { n , 1 , 100 } , { m , 1 , 100 } ] , 1 ] ;

t4 = Map [ f , t3 ] ;

p3 = ListPlot [ t4 , PlotRange → { { 0 , 100 } , { 0 , 100 } } , AspectRatio → 1 ,

PlotStyle → CMYKColor [ 0 , 1 , . 65 , . 34 ] ,

PlotLabel → “(c)” ] ;

With Mathematica, we can modify commands to investigate how changing parameters affect a given situation. In the following, we compute the sine of p if x=p/q. See Fig. 3.5 (d).

Clear [ f ]

f [ { a _ , b _ } ] := { a / b , Sin [ ( a / b // Numerator ) ] } ;

t5 = Flatten [ Table [ { n , m } , { n , 1 , 300 } , { m , 1 , 200 } ] , 1 ] ;

t6 = Map [ f , t5 ] ;

p4 = ListPlot [ t6 , PlotRange → { { 0 , 3 } , { 0 , 1 } } , AspectRatio → 1 ,

PlotStyle → CMYKColor [ 0 , . 65 , 1 , . 09 ] ,

PlotLabel → “(d)” ] ;

Show[GraphicsGrid[{{p1,p2},{p3,p4}}]]  □

3.2 Differential Calculus

3.2.1 Definition of the Derivative

Definition 3.2

The derivative of y=f(x) is

y ′ = f ′ ( x ) = d y d x = lim h → 0 ⁡ f ( x + h ) − f ( x ) h ,

(3.1)

provided the limit exists.

Assuming the derivative exists, as h approaches 0 the secants approach the tangent. Hence, if the limit exists the derivative gives us the slope of a function at that particular value of x.

The Limit command can be used along with Simplify to compute the derivative of a function using the definition of the derivative.

Example 3.8

Use the definition of the derivative to compute the derivative of (a) f(x)=x+1/x, (b) g(x)=1/x and (c) f(x)=sin⁡2x.

Solution

For (a), we first define f, compute the difference quotient, (f(x+h)−f(x))/h, simplify the difference quotient with Simplify, and use Limit to calculate the derivative.

f [ x _ ] = x + 1 / x ;

step1 = ( f [ x + h ] − f [ x ] ) / h

step2 = Simplify [ step1 ]

Limit [ step2 , h → 0 ]

h − 1 x + 1 h + x h

− 1 + h x + x 2 x ( h + x )

1 − 1 x 2

For (b), we use the same approach as in (a) but use Together rather than Simplify to reduce the complex fraction.

g [ x _ ] = 1 / Sqrt [ x ] ;

step1 = ( g [ x + h ] − g [ x ] ) / h

step2 = Together [ step1 ]

Limit [ step2 , h → 0 ]

− 1 x + 1 h + x h

x − h + x h x h + x

− 1 2 x 3 / 2

For (c), observe that Simplify instructs Mathematica to apply elementary trigonometric identities.

f [ x _ ] = Sin [ 2 x ] ;

step1 = ( f [ x + h ] − f [ x ] ) / h

step2 = Simplify [ step1 ]

Limit [ step2 , h → 0 ]

−Sin[2x]+Sin[2(h+x)]h

2 Cos [ h + 2 x ] Sin [ h ] h

2Cos[2x] □

If the derivative of y=f(x) exists at x=a, a geometric interpretation of f′(a) is that f′(a) is the slope of the line tangent to the graph of y=f(x) at the point (a,f(a)).

To motivate the definition of the derivative, many calculus texts choose a value of x, x=a, and then draw the graph of the secant line passing through the points (a,f(a)) and (a+h,f(a+h)) for “small” values of h to show that as h approaches 0, the secant line approaches the tangent line. An equation of the secant line passing through the points (a,f(a)) and (a+h,f(a+h)) is given by

y − f ( a ) = f ( a + h ) − f ( a ) ( a + h ) − a ( x − a ) or y = f ( a + h ) − f ( a ) h ( x − a ) + f ( a ) .

Example 3.9

If f(x)=x2−4x, graph f(x) together with the secant line containing (1,f(1)) and (1+h,f(1+h)) for various values of h.

Solution

We begin by considering a particular h value. We choose h=0.4. We then define f(x)=x2−4x. In p1, we graph f(x) in black on the interval [−1,5], in p2 we place a blue point at (1,f(1)) and a green point at (1.4,f(1.4), in p3 we graph the tangent to y=f(x) at (1,f(1)) in red, in p4 we graph the secant containing (1,f(1)) and (1.4,f(1.4) in purple, and finally we show all four graphics together with Show in Fig. 3.6.

Figure 3.6 Plots of y = x² − 4x, the tangent at (1,f(1)), and the secant containing (1,f(1)) and (1 + h,f(1 + h)) if h = 0.4 (Georgia Tech colors).

Remember that when a semi-colon is placed at the end of a command, the resulting output is not displayed. The names of the colors that Mathematica knows are listed in the ColorSchemes palette followed by “Known” and then “System.”

f [ x _ ] = x ∧ 2 − 4 x ;

p1 = Plot [ f [ x ] , { x , − 1 , 5 } , PlotStyle → CMYKColor [ . 08 , . 3 , 1 , 0 ] ] ;

p2 = Graphics [ { PointSize [ . 03 ] , CMYKColor [ 1 , . 88 , . 39 , . 42 ] , Point [ { 1 , f [ 1 ] } ] ,

CMYKColor [ . 04 , . 14 , . 67 , 0 ] , Point [ { 1 + . 4 , f [ 1 + . 4 ] } ] } ] ;

p3 = Plot [ f ′ [ 1 ] ( x − 1 ) + f [ 1 ] , { x , − 1 , 5 } , PlotStyle → CMYKColor [ 1 , . 88 , . 39 , . 42 ] ] ;

p4 = Plot [ ( f [ 1 + . 4 ] − f [ 1 ] ) / . 4 ( x − 1 ) + f [ 1 ] , { x , − 1 , 5 } ,

PlotStyle → CMYKColor [ . 04 , . 14 , . 67 , 0 ] ] ;

Show [ p1 , p2 , p3 , p4 , PlotRange → { { − 1 , 5 } , { − 6 , 6 } } , AspectRatio → 1 ]

We now generalize the previous set of commands for arbitrary h≠0 values. g(h) shows plots of y=x2−4x, the tangent at (1,f(1)), and the secant containing (1,f(1)) and (1+h,f(1+h)).

Clear [ f , g ] ;

f [ x _ ] = x ∧ 2 − 4 x ;

g [ h _ ] := Module [ { p1 , p2 , p3 , p4 } ,

p1 = Plot [ f [ x ] , { x , − 1 , 5 } , PlotStyle → CMYKColor [ . 08 , . 3 , 1 , 0 ] ] ;

p2 = Graphics [ { PointSize [ . 03 ] , CMYKColor [ 1 , . 88 , . 39 , . 42 ] , Point [ { 1 , f [ 1 ] } ] ,

CMYKColor [ . 04 , . 14 , . 67 , 0 ] , Point [ { 1 + h , f [ 1 + h ] } ] } ] ;

p3 = Plot [ f ′ [ 1 ] ( x − 1 ) + f [ 1 ] , { x , − 1 , 5 } , PlotStyle → CMYKColor [ 1 , . 88 , . 39 , . 42 ] ] ;

p4 = Plot [ ( f [ 1 + h ] − f [ 1 ] ) / h ( x − 1 ) + f [ 1 ] , { x , − 1 , 5 } ,

PlotStyle - > CMYKColor [ . 04 , . 14 , . 67 , 0 ] ] ;

Show [ p1 , p2 , p3 , p4 , PlotRange → { { − 1 , 5 } , { − 6 , 6 } } , AspectRatio → 1 ] ]

Table[f[x],{x,start,stop,stepsize}] creates a table of f(x) values beginning with start and ending with stop using increments of stepsize. Given a table, Partition[table,n] partitions the table into n element subgroups. Thus if a table, t1 has 9 elements, Partition[t1, 3] creates a 3×3 grid; three sets of three elements each.

Using Table followed by GraphicsGrid, we can create an table of graphics for various values of h like that shown in Fig. 3.7. With Table, the dimensions of the grid displayed on your computer are based on the size of the active Mathematica window. To control the dimensions of the grid, we use GraphicsGrid together with Partition and Show.

Figure 3.7 Plots of y = x² − 4x, the tangent at (1,f(1)), and the secant containing (1,f(1)) and (1 + h,f(1 + h)) for various values of h.

t1 = Table [ g [ k ] , { k , 1 , . 0001 , − ( 1 − . 0001 ) / 8 } ] ;

Show [ GraphicsGrid [ Partition [ t1 , 3 ] ] ]

Animate works in the same way as Table and is very similar to Manipulate. Entering

Animate [ g [ k ] , { k , 1 , . 0001 , − ( 1 − . 0001 ) / 99 } ]

generates an animation of g(k) and 100 equally spaced values of k starting with k=1 and ending with k=0.0001. To animate the result, use the toolbar in the animate graphic. To vary g continuously and not use a specific stepsize enter Animate[g[k],{k,1,.0001}].

After animating the selection, you can control the animation (speed, direction, pause, and so on) with the buttons at the top of the animation window.

With Mathematica 11, you can use Manipulate to help generate animations and images that you can adjust based on changing parameter values.

To illustrate how to do so, we begin by redefining f and then defining m(a,h). Given a and h values, m(a,h) plots f(x) for −10⩽x⩽10 (p1), plots a blue point at (a,f(a)) and a green point at (a+h,f(a+h)) (p2), plots f′(a)(x−a)+f(a) (the tangent to the graph of f(x) at (a,f(a))) for −1⩽x⩽5 in red (p3), the secant containing (a,f(a)) and (a+h,f(a+h)) for −10⩽x⩽10 in purple (p4) and finally displays all four graphics together with Show. Using PlotRange, we indicate that horizontal axis displays x values between −10 and 10, the vertical axis displays y values between −10 and 10; AspectRatio->1 means that the ratio of the lengths of the x to y axes is 1. Thus, the plot scaling is correct. Note that when we use Module to define m, p1, p2, p3, and p4 are local to the function m. This means that if you have such objects defined elsewhere in your Mathematica notebook, those objects are not affected when you compute m.

Clear [ m , f ] ;

f [ x _ ] = x ∧ 2 − 4 x ;

m [ a _ , h _ ] := Module [ { p1 , p2 , p3 , p4 } ,

p1 = Plot [ f [ x ] , { x , − 10 , 10 } , PlotStyle → CMYKColor [ . 08 , . 3 , 1 , 0 ] ] ;

p2 = Graphics [ { PointSize [ . 03 ] , CMYKColor [ 1 , . 88 , . 39 , . 42 ] , Point [ { a , f [ a ] } ] ,

CMYKColor [ . 04 , . 14 , . 67 , 0 ] , Point [ { a + h , f [ a + h ] } ] } ] ;

p3 = Plot [ f ′ [ a ] ( x − a ) + f [ a ] , { x , − 1 , 5 } , PlotStyle → CMYKColor [ 1 , . 88 , . 39 , . 42 ] ] ;

p4 = Plot [ ( f [ a + h ] − f [ a ] ) / h ( x − a ) + f [ a ] , { x , − 10 , 10 } ,

PlotStyle → CMYKColor [ . 04 , . 14 , . 67 , 0 ] ] ;

Show [ p1 , p2 , p3 , p4 , PlotRange → { { − 10 , 10 } , { − 10 , 10 } } , AspectRatio → 1 ] ]

Now we use Manipulate to create a “mini” program. The sliders (centered at a=0 and h=.5 with range from −10 to 10 an −1 to 1, respectively) allow you to see how changing a and h affects the plot. See Fig. 3.8.

Figure 3.8 With Manipulate, we can perform animations and see how a function changes depending on parameter values. When you press the + button, the bar expands so that you can control the Manipulate object.

Manipulate [ m [ a , h ] , { { a , 0 } , − 10 , 10 } , { { h , . 5 } , − 1 , 1 } ]

Fig. 3.8 illustrates the special case when f(x)=x2−4x. To illustrate the same concept using a “standard” set of functions (polynomials, rational, root, and trig), we first define the functions

quad [ x _ ] = ( x + 2 ) ∧ 2 − 2 ;

cubic [ x _ ] = − 1 / 10 x ( x ∧ 2 − 25 ) ;

rational [ x _ ] = 50 / ( ( x + 5 ) ( x − 5 ) ) ;

root [ x _ ] = 3 Sqrt [ x + 5 ] ;

sin ⁡ [ x _ ] = 5 Sin [ x ] ;

and then we adjust m by defining a few of these “standard” and then defining the function mmore, which performs the same actions as m but does so for the function selected. We then use Manipulate to create an object that shows the secant (in yellow), the tangent (in blue) for the selected function, a value, and h value. See Fig. 3.9.

Figure 3.9 With this Manipulate object, we see how various functions, a values, and h values affect the secant to y = f(x) passing through (a,f(a)) and (a + h,f(a + h)) and the tangent to y = f(x) at (a,f(a)). (For interpretation of the colors in this figure, the reader is referred to the web version of this chapter.)

Clear [ mmore ] ;

mmore [ f _ , a _ , h _ ] := Module [ { p1 , p2 , p3 , p4 } ,

p1 = Plot [ f [ x ] , { x , − 10 , 10 } , PlotStyle → CMYKColor [ . 08 , . 3 , 1 , 0 ] ] ;

p2 = Graphics [ { PointSize [ . 03 ] , CMYKColor [ 1 , . 88 , . 39 , . 42 ] , Point [ { a , f [ a ] } ] ,

CMYKColor [ . 04 , . 14 , . 67 , 0 ] , Point [ { a + h , f [ a + h ] } ] } ] ;

p3 = Plot [ f ′ [ a ] ( x − a ) + f [ a ] , { x , − 10 , 10 } , PlotStyle → CMYKColor [ 1 , . 88 , . 39 , . 42 ] ] ;

p4 = Plot [ ( f [ a + h ] − f [ a ] ) / h ( x − a ) + f [ a ] , { x , − 10 , 10 } ,

PlotStyle → CMYKColor [ . 04 , . 14 , . 67 , 0 ] ] ;

Show [ p1 , p2 , p3 , p4 , PlotRange → { { − 10 , 10 } , { − 10 , 10 } } , AspectRatio → 1 ] ]

Manipulate [ mmore [ f , a , h ] , { { f , quad } , { quad , cubic , rational , root , sin ⁡ } } ,

{{a,0},−10,10},{{h,1},−2,2}]  □

3.2.2 Calculating Derivatives

The functions D and ' are used to differentiate functions. Assuming that y=f(x) is differentiable,

1.  D[f[x],x] computes and returns f′(x)=df/dx,

2.  f'[x] computes and returns f′(x)=df/dx,

3.  f''[x] computes and returns f(2)(x)=d2f(x)/dx2, and

4.  D[f[x],{x,n}] computes and returns f(n)(x)=dnf(x)/dxn.

5.  You can use the button located on the Basic Math Assistant and Basic Math Input palettes to create templates to compute derivatives.

Fig. 3.10 illustrates various ways of computing derivatives using the ′ symbol, D, and the ∂ symbol.

Figure 3.10 You can use ′, D, and ∂ to compute derivatives of functions.

Mathematica knows the numerous differentiation rules, including the product, quotient, and chain rules. Thus, entering

Clear [ f , g ]

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

g [ x ] f ′ [ x ] + f [ x ] g ′ [ x ]

shows us that ddx(f(x)⋅g(x))=f′(x)g(x)+f(x)g′(x); entering

Together [ D [ f [ x ] / g [ x ] , x ] ]

g [ x ] f ′ [ x ] − f [ x ] g ′ [ x ] g [ x ] 2

shows us that

Remark 3.4

Throughout the text, input is in Bold and output is Not; output follows input ddx(f(x)/g(x))=(f′(x)g(x)−f(x)g′(x))/(g(x))2; and entering

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

shows us that ddx(f(g(x)))=f′(g(x))g′(x).

Example 3.10

Compute the first and second derivatives of (a) y=x4+43x3−3x2, (b) f(x)=4x5−52x4−10x3, (c) y=e2x+e−2x, and (d) y=(1+1/x)x.

Solution

For (a), we use D.

D[x∧4+4/3x∧3−3x∧2,x] −6x+4x2+4x3

For (b), we first define f and then use ′ together with Factor to calculate and factor f′(x) and f″(x).

f [ x _ ] = 4 x ∧ 5 − 5 / 2 x ∧ 4 − 10 x ∧ 3 ;

Factor[f′[x]] 10x2(1+x)(−3+2x) Factor[f″[x]] 10x(−6−3x+8x2)

For (c), we use Simplify together with D to calculate and simplify y′ and y″.

D[Sqrt[Exp[2x]+Exp[−2x]],x]

− 2 e − 2 x + 2 e 2 x 2 e − 2 x + e 2 x

D [ Sqrt [ Exp [ 2 x ] + Exp [ − 2 x ] ] , { x , 2 } ] // Simplify

e − 2 x + e 2 x ( 1 + 6 e 4 x + e 8 x ) ( 1 + e 4 x ) 2

By hand, (d) would require logarithmic differentiation. The second derivative would be particularly difficult to compute by hand. Mathematica quickly computes and simplifies each derivative.

Simplify [ D [ ( 1 + 1 / x ) ∧ x , x ] ]

( 1 + 1 x ) x ( − 1 + ( 1 + x ) Log [ 1 + 1 x ] ) 1 + x

Simplify [ D [ ( 1 + 1 / x ) ∧ x , { x , 2 } ] ]

(1+1x)x(−1+x−2x(1+x)Log[1+1x]+x(1+x)2Log[1+1x]2)x(1+x)2  □

The command Map[f,list] applies the function f to each element of the list list. Thus, if you are computing the derivatives of a large number of functions, you can use Map together with D.

Map and operations on lists are discussed in more detail in Chapter 4.

A built-in Mathematica function is threadable if f[list] returns the same result as Map[f,list]. Many familiar functions like D and Integrate are threadable.

Example 3.11

Compute the first and second derivatives of sin⁡x, cos⁡x, tan⁡x, sin−1⁡x, cos−1⁡x, and tan−1⁡x.

Solution

Notice that lists are contained in braces. Thus, entering

Map [ D [ # , x ] & , { Sin [ x ] , Cos [ x ] , Tan [ x ] ,

ArcSin [ x ] , ArcCos [ x ] , ArcTan [ x ] } ]

{ Cos [ x ] , − Sin [ x ] , Sec [ x ] 2 , 1 1 − x 2 , − 1 1 − x 2 , 1 1 + x 2 }

computes the first derivative of the three trigonometric functions and their inverses. In this case, we have applied a pure function to the list of trigonometric functions and their inverses. Given an argument #, D[#,x]& computes the derivative of # with respect to x. The & symbol is used to mark the end of a pure function. Similarly, entering

Map [ D [ # , { x , 2 } ] & , { Sin [ x ] , Cos [ x ] , Tan [ x ] ,

ArcSin [ x ] , ArcCos [ x ] , ArcTan [ x ] } ]

{−Sin[x],−Cos[x],2Sec[x]2Tan[x],x(1−x2)3/2,−x(1−x2)3/2,−2x(1+x2)2}

computes the second derivative of the three trigonometric functions and their inverses. Because D is threadable, the same results are obtained with the following commands.

D [ { Sin [ x ] , Cos [ x ] , Tan [ x ] ,

ArcSin [ x ] , ArcCos [ x ] , ArcTan [ x ] } , x ]

{ Cos [ x ] , − Sin [ x ] , Sec [ x ] 2 , 1 1 − x 2 , − 1 1 − x 2 , 1 1 + x 2 }

D [ { Sin [ x ] , Cos [ x ] , Tan [ x ] ,

ArcSin [ x ] , ArcCos [ x ] , ArcTan [ x ] } , { x , 2 } ]

{−Sin[x],−Cos[x],2Sec[x]2Tan[x],x(1−x2)3/2,−x(1−x2)3/2,−2x(1+x2)2} □

With DynamicModule, we create a simple dynamic that lets you compute the first and second derivatives of basic functions and plot them on a standard viewing window, [−5,5]×[−5,5]. The layout of Fig. 3.11 is primarily determined by Panel, Column, and Grid. The default function is y=x2. To compute and graph the first and second derivatives of a different function, simply type over x2 with the desired function.

Figure 3.11 Seeing the relationship between the first and second derivative of a function and the original function (Mercer University colors).

Panel [ DynamicModule [ { f = x ∧ 2 } ,

Column [ { InputField [ Dynamic [ f ] ] , Grid [ { { “First Derivative” ,

Panel [ Dynamic [ D [ f , x ] // Simplify ] ] } ,

{ “Second Derivative” , Panel [ Dynamic [ D [ f , { x , 2 } ] // Simplify ] ] } } ] ,

Dynamic [ Plot [ Evaluate [ Tooltip [ { f , D [ f , x ] , D [ f , { x , 2 } ] } ] ] ,

{ x , − 5 , 5 } , PlotRange → { − 5 , 5 } ,

PlotStyle → { { CMYKColor [ 0 , . 65 , 1 , . 04 ] } , { Black } , { CMYKColor [ . 01 , . 05 , . 23 , . 03 ] } } ,

AspectRatio → Automatic ] ] } ] ] , ImageSize → { 300 , 300 } ]

3.2.3 Implicit Differentiation

If an equation contains two variables, x and y, implicit differentiation can be carried out by explicitly declaring y to be a function of x, y=y(x), and using D or by using the Dt command.

Example 3.12

Find y′=dy/dx if (a) cos⁡(exy)=x and (b) ln⁡(x/y)+5xy=3y.

Solution

For (a) we illustrate the use of D. Notice that we are careful to specifically indicate that y=y(x). First we differentiate with respect to x

Clear [ x , y ]

s1 = D [ Cos [ Exp [ x y [ x ] ] ] − x , x ]

−1−exy[x]Sin[exy[x]](y[x]+xy′[x]) and then we solve the resulting equation for y′=dy/dx with Solve.

Solve [ s1 = = 0 , y ′ [ x ] ]

{ { y ′ [ x ] → − e − x y [ x ] ( Csc [ e x y [ x ] ] + e x y [ x ] y [ x ] ) x } }

For (b), we use Dt. When using Dt, we interpret Dt[x]=1 and Dt[y]=y′=dy/dx. Thus, entering

s2 = Dt [ Log [ x / y ] + 5 x y − 3 y ]

5 y Dt [ x ] − 3 Dt [ y ] + 5 x Dt [ y ] + y ( Dt [ x ] y − x Dt [ y ] y 2 ) x

s3 = s2 /. { Dt [ x ] → 1 , Dt [ y ] → dydx }

− 3 dydx + 5 dydx x + 5 y + ( − dydx x y 2 + 1 y ) y x

and solving for dydx with Solve

Solve [ s3 = = 0 , dydx ]

{ { dydx → − y ( 1 + 5 x y ) x ( − 1 − 3 y + 5 x y ) } }

shows us that if ln⁡(x/y)+5xy=3y, y′=dydx=−(1+5xy)y(5xy−3y−1)x.

To graph each equation, we use ContourPlot. Generally, given an equation of the form f(x,y)=g(x,y), the command

ContourPlot[f[x,y]==g[x,y],{x,a,b},{y,c,d}]

attempts to plot the graph of f(x,y)=g(x,y) on the rectangle [a,b]×[c,d]. Using Show together with GraphicsRow, we show the two graphs side-by-side in Fig. 3.12.

Figure 3.12 On the left, cos⁡(exy)=x for −2 ⩽ x ⩽ 2 and −4 ⩽ y ⩽ 4; on the right, ln⁡(x/y)+5xy=3y for .01 ⩽ x ⩽ 3 and .01 ⩽ y ⩽ 3 (MIT colors).

cp1 = ContourPlot [ Cos [ Exp [ x y ] ] == x , { x , − 2 , 2 } , { y , − 4 , 4 } , PlotPoints → 120 ,

Frame → False , Axes → Automatic , AxesOrigin → { 0 , 0 } , Contours → 40 ,

ContourStyle → CMYKColor [ . 24 , 1 , . 78 , . 17 ] , AxesLabel → { x , y } ,

PlotLabel → “(a)” ] ;

cp2 = ContourPlot [ Log [ x / y ] + 5 x y = = 3 y , { x , . 01 , 3 } , { y , . 01 , 3 } , PlotPoints → 120 ,

AxesLabel → { x , y } , PlotLabel → “(b)” ,

ContourStyle → CMYKColor [ . 48 , . 39 , . 39 , . 04 ] ,

Frame → False , Axes → Automatic , AxesOrigin → { 0 , 0 } ] ;

Show[GraphicsRow[{cp1,cp2}]] □

3.2.4 Tangent Lines

If f′(a) exists, f′(a) is interpreted to be the slope of the line tangent to the graph of y=f(x) at the point (a,f(a)). In this case, an equation of the tangent is given by

y − f ( a ) = f ′ ( a ) ( x − a ) or y = f ′ ( a ) ( x − a ) + f ( a ) .

Example 3.13

Find an equation of the line tangent to the graph of f(x)=sin⁡x1/3+cos1/3⁡x at the point with x-coordinate x=5π/3.

Solution

Recall that when computing odd roots of negative numbers, Mathematica returns the value with the largest imaginary part. However, for our purposes, we need the real-valued root. To obtain the real-valued root, use Surd: Surd[x,n] returns the real-valued root of x if n is odd. Because we will be graphing a function involving odd roots of negative numbers, we define f(x) using the Surd function and then compute f′(x).

f [ x _ ] = Sin [ Surd [ x , 3 ] ] + Surd [ Cos [ x ] , 3 ] ;

f ′ [ x ]

Cos[x3]3x32−Sin[x]3Cos[x]32

Then, the slope of the line tangent to the graph of f(x) at the point with x-coordinate x=5π/3 is

f ′ [ 5 Pi / 3 ]

1 2 1 / 3 3 + Cos [ ( 5 π 3 ) 1 / 3 ] 3 1 / 3 ( 5 π ) 2 / 3

f ′ [ 5 Pi / 3 ] // N

0.440013

while the y-coordinate of the point is

f [ 5 Pi / 3 ]

1 2 1 / 3 + Sin [ ( 5 π 3 ) 1 / 3 ]

f [ 5 Pi / 3 ] // N

1.78001

Thus, an equation of the line tangent to the graph of f(x) at the point with x-coordinate x=5π/3 is

y − ( 1 2 3 + sin ⁡ 5 π / 3 3 ) = ( cos ⁡ 5 π / 3 3 3 3 25 π 2 3 + 1 2 3 3 ) ( x − 5 π 3 ) ,

as shown in Fig. 3.13.

Figure 3.13 f(x)=sin⁡x1/3+cos1/3⁡x together with its tangent at the point (5π/3,f(5π/3)) (Stanford colors).

p1 = Plot [ f [ x ] , { x , 0 , 4 Pi } , PlotStyle → CMYKColor [ 0 , 1 , . 65 , . 34 ] ] ;

p2 = ListPlot [ { { 5 Pi / 3 , f [ 5 Pi / 3 ] } // N } , PlotStyle → { PointSize [ . 03 ] ,

CMYKColor [ . 48 , . 36 , . 24 , . 66 ] } ] ;

p3 = Plot [ f ′ [ 5 Pi / 3 ] ( x − 5 Pi / 3 ) + f [ 5 Pi / 3 ] , { x , 0 , 4 Pi } ,

PlotStyle->CMYKColor[.48,.36,.24,.66]];

Show [ p1 , p2 , p3 , AspectRatio - > Automatic , PlotRange → All ,

AxesLabel → { x , y } ]

We use Manipulate to plot different tangents or animate the tangents as shown in Fig. 3.14.

Figure 3.14 Using Manipulate to plot f(x) together with various tangents.

Manipulate [

p1 = Plot [ f [ x ] , { x , 0 , 4 Pi } , PlotStyle → CMYKColor [ 0 , 1 , . 65 , . 34 ] ] ;

p2 = ListPlot [ { { a , f [ a ] } // N } , PlotStyle → { PointSize [ . 03 ] ,

CMYKColor [ . 48 , . 36 , . 24 , . 66 ] } ] ;

p3 = Plot [ f ′ [ a ] ( x − a ) + f [ a ] , { x , 0 , 4 Pi } , PlotStyle - > CMYKColor [ . 48 , . 36 , . 24 , . 66 ] ] ;

Show [ p1 , p2 , p3 , AspectRatio - > Automatic , PlotRange → { − 1 , 5 } ,

AxesLabel→{x,y}],{{a,5},.1,4Pi−.1}]  □

Remark 3.5

Sometimes using Surd may be cumbersome. Therefore, it is worth noting that if you desire a real root to xn/m and there is a real root, x^(n/m)=Abs[x]^(n/m) if n is even and x^(n/m)=Sign[x]Abs[x]^(n/m) if n is odd. Of course, remember that if m is even and x is negative, xn/m is not a real number. Note that Sign[x] returns 1 if x is positive and −1 if x is negative.

Tangent Lines of Implicit Functions

Example 3.14

Find equations of the tangent line and normal line to the graph of x2y−y3=8 at the point (−3,1). Find and simplify y″=d2y/dx2.

If the line ℓ has slope m1 the line perpendicular (or, normal) to ℓ has slope −1/m1.

Solution

We evaluate y′=dy/dx if x=−3 and y=1 to determine the slope of the tangent line at the point (−3,1). Note that we cannot (easily) solve x2y−y3=8 for y so we use implicit differentiation to find y′=dy/dx:

d d x ( x 2 y − y 3 ) = d d x ( 8 )

2 x y + x 2 y ′ − 3 y 2 y ′ = 0 y ′ = − 2 x y x 2 − 3 y 2 .

By the product and chain rules, ddx(x2y)=ddx(x2)y+x2ddx(y)=2x⋅y+x2⋅dydx=2xy+x2y′.

eq = x ∧ 2 y − y ∧ 3 = = 8

x 2 y − y 3 = = 8

s1 = Dt [ eq ]

2 x y Dt [ x ] + x 2 Dt [ y ] − 3 y 2 Dt [ y ] = = 0

s2 = s1 /. Dt [ x ] → 1

2 x y + x 2 Dt [ y ] − 3 y 2 Dt [ y ] = = 0

s3 = Solve [ s2 , Dt [ y ] ]

{ { Dt [ y ] → − 2 x y x 2 − 3 y 2 } }

Notice that s3 is a list. The formula for y′=dy/dx is the second part of the first part of the first part of s3 and extracted from s3 with

Lists are discussed in more detail in Chapter 4.

s3 [ [ 1 , 1 , 2 ] ]

− 2 x y x 2 − 3 y 2

We then use ReplaceAll (/.) to find that the slope of the tangent at (−3,1) is

s3 [ [ 1 , 1 , 2 ] ] /. { x → − 3 , y → 1 }

1

The slope of the normal is −1/1=−1. Equations of the tangent and normal are given by

y − 1 = 1 ( x + 3 ) and y − 1 = − 1 ( x + 3 ) ,

respectively. See Fig. 3.15.

Figure 3.15 Graphs of x²y − y³ = 8 and the tangent and normal at (−3,1) (University of California Berkeley colors).

cp1 = ContourPlot [ x ∧ 2 y − y ∧ 3 − 8 , { x , − 5 , 5 } , { y , − 5 , 5 } , Contours → { 0 } ,

ContourShading → False , PlotPoints → 200 ,

ContourStyle → { Thickness [ . 01 ] , CMYKColor [ 1 , . 45 , 0 , . 66 ] } ] ;

p1 = ListPlot [ { { − 3 , 1 } } , PlotStyle → { Black , PointSize [ . 03 ] } ] ;

p2 = Plot [ { ( x + 3 ) + 1 , − ( x + 3 ) + 1 } , { x , − 5 , 5 } ,

PlotStyle → { { CMYKColor [ 0 , . 3 , . 93 , 0 ] , Thickness [ . 01 ] } } ] ;

Show [ cp1 , p1 , p2 , Frame → False , Axes → Automatic , AxesOrigin → { 0 , 0 } ,

AspectRatio → Automatic , DisplayFunction → $ DisplayFunction ]

To find y″=d2y/dx2, we proceed as follows.

s4 = Dt [ s3 [ [ 1 , 1 , 2 ] ] ] // Simplify

− 2 ( x 2 + 3 y 2 ) ( − y Dt [ x ] + x Dt [ y ] ) ( x 2 − 3 y 2 ) 2

s5 = s4 /. Dt [ x ] → 1 /. s3 [ [ 1 ] ] // Simplify

6 y ( x 2 − y 2 ) ( x 2 + 3 y 2 ) ( x 2 − 3 y 2 ) 3

The result means that

y ″ = d 2 y d x 2 = 6 ( x 2 y − y 3 ) ( x 2 + 3 y 2 ) ( x 2 − 3 y 2 ) 3 .

Because x2y−y3=8, the second derivative is further simplified to

y ″ = d 2 y d x 2 = 48 ( x 2 + 3 y 2 ) ( x 2 − 3 y 2 ) 3 .

 □

Parametric Equations and Polar Coordinates

For the parametric equations {x=f(t),y=g(t)}, t∈I,

y ′ = d y d x = d y / d t d x / d t = g ′ ( t ) f ′ ( t )

and

y ″ = d 2 y d x 2 = d d x d y d x = d / d t ( d y / d x ) d x / d t .

If {x=f(t),y=g(t)} has a tangent line at the point (f(a),g(a)), parametric equations of the tangent are given by

x = f ( a ) + t f ′ ( a ) and y = g ( a ) + t g ′ ( a ) .

(3.2)

If f′(a), g′(a)≠0, we can eliminate the parameter from (3.2)

x − f ( a ) f ′ ( a ) = y − g ( a ) g ′ ( a ) y − g ( a ) = g ′ ( a ) f ′ ( a ) ( x − f ( a ) )

and obtain an equation of the tangent line in point-slope form.

l = Solve [ x [ a ] + t x ′ [ a ] = = cx , t ]

r = Solve [ y [ a ] + t y ′ [ a ] = = cy , t ]

{ { t → cx − x [ a ] x ′ [ a ] } }

{ { t → cy − y [ a ] y ′ [ a ] } }

Example 3.15

The Cycloid

The cycloid has parametric equations

x = t − sin ⁡ t and y = 1 − cos ⁡ t .

Graph the cycloid together with the line tangent to the graph of the cycloid at the point (x(a),y(a)) for various values of a between −2π and 4π.

Solution

After defining x and y we use ' to compute dy/dt and dx/dt. We then compute dy/dx=(dy/dt)/(dx/dt) and d2y/dx2.

x [ t _ ] = t − Sin [ t ] ;

y [ t _ ] = 1 − Cos [ t ] ;

dx = x ′ [ t ]

dy = y ′ [ t ]

dydx = dy / dx

1 − Cos [ t ]

Sin [ t ]

Sin [ t ] 1 − Cos [ t ]

dypdt=Simplify[D[dydx,t]]

1 − 1 + Cos [ t ]

secondderiv = Simplify [ dypdt / dx ]

− 1 ( − 1 + Cos [ t ] ) 2

We then use ParametricPlot to graph the cycloid for −2π⩽t⩽4π, naming the resulting graph p1.

p1 = ParametricPlot [ { x [ t ] , y [ t ] } , { t , − 2 Pi , 4 Pi } ,

PlotStyle → { { CMYKColor [ 0 , . 68 , . 98 , 0 ] , Thickness [ . 015 ] } } ] ;

Next, we use Table to define toplot to be 40 tangent lines (3.2) using equally spaced values of a between −2π and 4π. We then graph each line toplot and name the resulting graph p2. Finally, we show p1 and p2 together with the Show function. The resulting plot is shown to scale because the lengths of the x and y-axes are equal and we include the option AspectRatio->1. In the graphs, notice that on intervals for which dy/dx is defined, dy/dx is a decreasing function and, consequently, d2y/dx2<0. (See Fig. 3.16.)

Figure 3.16 The cycloid with various tangents (California Institute of Technology colors).

toplot = Table [ { x [ a ] + t x ′ [ a ] , y [ a ] + t y ′ [ a ] } , { a , − 2 Pi , 4 Pi , 6 Pi / 39 } ] ;

p2 = ParametricPlot [ Evaluate [ toplot ] , { t , − 2 , 2 } , PlotStyle → CMYKColor [ . 35 , . 28 , . 35 , 0 ] ] ;

Show [ p1 , p2 , AspectRatio → 1 , PlotRange → { − 3 Pi , 3 Pi } ,

AxesLabel → { x , y } ]

With Manipulate, you can animate the tangents. (See Fig. 3.17.)

Figure 3.17 With Manipulate or Animate you can animate the tangents.

Manipulate [ x [ t _ ] = t − Sin [ t ] ;

y [ t _ ] = 1 − Cos [ t ] ;

y [ t _ ] = Module [ { p1 , p2 } , p1 = ParametricPlot [ { x [ t ] , y [ t ] } , { t , − 2 Pi , 4 Pi } ,

PlotStyle → { { CMYKColor [ 0 , . 68 , . 98 , 0 ] , Thickness [ . 01 ] } } ] ;

p2 = ParametricPlot [ { x [ a ] + t x ′ [ a ] , y [ a ] + t y ′ [ a ] } , { t , − 5 , 5 } ,

PlotStyle → { { CMYKColor [ . 35 , . 28 , . 35 , 0 ] , Thickness [ . 0075 ] } } ] ;

Show [ p1 , p2 , AspectRatio → 1 , PlotRange → { { − 2 Pi , 4 Pi } , { − 3 Pi , 3 Pi } } ] ] ,

{{a,1},−2Pi,4Pi}]  □

Example 3.16

Orthogonal Curves

Two lines L1 and L2 with slopes m1 and m2, respectively, are orthogonal if their slopes are negative reciprocals: m1=−1/m2.

Extended to curves, we say that the curves C1 and C2 are orthogonal at a point of intersection if their respective tangent lines to the curves at that point are orthogonal.

Show that the family of curves with equation x2+2xy−y2=C is orthogonal to the family of curves with equation y2+2xy−x2=C.

Solution

We begin by defining eq1 and eq2 to be equations x2+2xy−y2=C and y2+2xy−x2=C, respectively. Then, use Dt to differentiate and Solve to find y′=dy/dx.

eq1 = x 2 + 2 x y − y 2 == c ;

eq2 = y 2 + 2 x y − x 2 == c ;

Simplify [ Solve [ Dt [ eq1 , x ] , Dt [ y , x ] ] /. Dt [ c , x ] → 0 ]

{{Dt[y,x]→−x+yx−y}}

Simplify [ Solve [ Dt [ eq2 , x ] , Dt [ y , x ] ] /. Dt [ c , x ] → 0 ]

{ { Dt [ y , x ] → x − y x + y } }

Because the derivatives are negative reciprocals, we conclude that the curves are orthogonal. We confirm this graphically by graphing several members of each family with ContourPlot and showing the results together. (See Fig. 3.18.)

Figure 3.18 x² + 2xy − y² = C and y² + 2xy − x² = C for various values of C (University of Illinois colors).

cp1 = ContourPlot [ x 2 + 2 x y − y 2 , { x , − 5 , 5 } , { y , − 5 , 5 } , ContourShading → False ,

ContourStyle → { { Thickness [ . 01 ] , CMYKColor [ 1 , . 9 , . 1 , . 77 ] } } ] ;

cp2 = ContourPlot [ y 2 + 2 x y − x 2 , { x , − 5 , 5 } , { y , − 5 , 5 } , ContourShading → False ,

ContourStyle → { { Thickness [ . 01 ] , CMYKColor [ 0 , . 76 , 1.0 ] } } ] ;

Show[cp1,cp2,Frame→False,Axes→Automatic,AxesOrigin→{0,0}]  □

Theorem 3.1

The Mean-Value Theorem for Derivatives

If y=f(x) is continuous on [a,b] and differentiable on (a,b) then there is at least one value of c between a and b for which

f ′ ( c ) = f ( b ) − f ( a ) b − a or, equivalently, f ( b ) − f ( a ) = f ′ ( c ) ( b − a ) .

(3.3)

Example 3.17

Find all number(s) c that satisfy the conclusion of the Mean-Value Theorem for f(x)=x2−3x on the interval [0,7/2].

Solution

By the power rule, f′(x)=2x−3. The slope of the secant containing (0,f(0)) and (7/2,f(7/2)) is

f ( 7 / 2 ) − f ( 0 ) 7 / 2 − 0 = 1 2 .

Solving 2x−3=1/2 for x gives us x=7/4.

f [ x _ ] = x ∧ 2 − 3 x

− 3 x + x 2

Solve [ f ′ [ x ] = = 0 , x ]

{ { x → 3 2 } }

Solve [ f ′ [ x ] = = ( f [ 7 / 2 ] − f [ 0 ] ) / ( 7 / 2 − 0 ) ]

{ { x → 7 4 } }

x=7/4 satisfies the conclusion of the Mean-Value Theorem for f(x)=x2−3x on the interval [0,7/2], as shown in Fig. 3.19.

Figure 3.19 Graphs of f(x)=x² − 3x, the secant containing (0,f(0)) and (7/2,f(7/2)), and the tangent at (7/4,f(7/4)) (University of Michigan colors).

p1 = Plot [ f [ x ] , { x , − 1 , 4 } , PlotStyle → CMYKColor [ 1 , . 6 , 0 , . 6 ] ] ;

p2 = Plot [ f [ x ] , { x , 0 , 7 / 2 } , PlotStyle → { { Thickness [ . 015 ] , CMYKColor [ 0 , . 18 , 1 , 0 ] } } ] ;

p3 = ListPlot [ { { 0 , f [ 0 ] } , { 7 / 4 , f [ 7 / 4 ] } , { 7 / 2 , f [ 7 / 2 ] } } ,

PlotStyle → { { CMYKColor [ . 2 , . 14 , . 12 , . 4 ] , PointSize [ . 025 ] } } ] ;

p4 = Plot [ { f ′ [ 7 / 4 ] ( x − 7 / 4 ) + f [ 7 / 4 ] , ( f [ 7 / 2 ] − f [ 0 ] ) / ( 7 / 2 − 0 ) x + f [ 0 ] } ,

{ x , − 2 , 4 } , PlotStyle → { { Dashing [ { . 01 } ] , CMYKColor [ . 06 , . 14 , . 38 , . 08 ] } ,

{ Dashing [ { . 01 } ] , CMYKColor [ . 21 , . 15 , . 54 , . 31 ] } } ] ;

Show [ p1 , p2 , p3 , p4 , DisplayFunction → $ DisplayFunction , AspectRatio → Automatic ,

PlotRange→{−3,3}]  □

3.2.5 The First Derivative Test and Second Derivative Test

Example 3.15 illustrates the following properties of the first and second derivative.

Theorem 3.2

Let y=f(x) be continuous on [a,b] and differentiable on (a,b).

1.  If f′(x)=0 for all x in (a,b), then f(x) is constant on [a,b].

2.  If f′(x)>0 for all x in (a,b), then f(x) is increasing on [a,b].

3.  If f′(x)<0 for all x in (a,b), then f(x) is decreasing on [a,b].

For the second derivative, we have the following theorem.

Theorem 3.3

Let y=f(x) have a second derivative on (a,b).

1.  If f″(x)>0 for all x in (a,b), then the graph of f(x) is concave up on (a,b).

2.  If f″(x)<0 for all x in (a,b), then the graph of f(x) is concave down on (a,b).

The critical points correspond to those points on the graph of y=f(x) where the tangent line is horizontal or vertical; the number x=a is a critical number if f′(a)=0 or f′(x) does not exist if x=a. The inflection points correspond to those points on the graph of y=f(x) where the graph of y=f(x) is neither concave up nor concave down. Theorems 3.2 and 3.3 help establish the first derivative test and second derivative test.

Theorem 3.4

First Derivative Test

Let x=a be a critical number of a function y=f(x) continuous on an open interval I containing x=a. If f(x) is differentiable on I, except possibly at x=a, f(a) can be classified as follows.

1.  If f′(x) makes a simple change in sign from positive to negative at x=a, then f(a) is a relative maximum.

2.  If f′(x) makes a simple change in sign from negative to positive at x=a, then f(a) is a relative minimum.

Theorem 3.5

Second Derivative Test

Let x=a be a critical number of a function y=f(x) and suppose that f″(x) exists on an open interval containing x=a.

1.  If f″(a)<0, then f(a) is a relative maximum.

2.  If f″(a)>0, then f(a) is a relative minimum.

Example 3.18

Graph f(x)=3x5−5x3.

Solution

We begin by defining f(x) and then computing and factoring f′(x) and f″(x).

f [ x _ ] = 3 x ∧ 5 − 5 x ∧ 3 ;

d1 = Factor [ f ′ [ x ] ]

d2 = Factor [ f ” [ x ] ]

15 ( − 1 + x ) x 2 ( 1 + x )

30 x ( − 1 + 2 x 2 )

By inspection, we see that the critical numbers are x=0, 1, and −1 while f″(x)=0 if x=0, 1/2, or −1/2. Of course, these values can also be found with Solve as done next in cns and ins, respectively.

cns = Solve [ d1 = = 0 ]

ins = Solve [ d2 = = 0 ]

{ { x → − 1 } , { x → 0 } , { x → 0 } , { x → 1 } }

{ { x → 0 } , { x → − 1 2 } , { x → 1 2 } }

We find the critical and inflection points by using /. (Replace All) to compute f(x) for each value of x in cns and ins, respectively. The result means that the critical points are (0,0), (1,−2) and (−1,2); the inflection points are (0,0), (1/2,−72/8), and (−1/2,72/8). We also see that f″(0)=0 so Theorem 3.5 cannot be used to classify f(0). On the other hand, f″(1)=30>0 and f″(−1)=−30<0 so by Theorem 3.5, f(1)=−2 is a relative minimum and f(−1)=2 is a relative maximum.

cps = { x , f [ x ] } /. cns

{ { − 1 , 2 } , { 0 , 0 } , { 0 , 0 } , { 1 , − 2 } }

f ” [ x ] /. cns

{ − 30 , 0 , 0 , 30 }

ips = { x , f [ x ] } /. ins

{ { 0 , 0 } , { − 1 2 , 7 4 2 } , { 1 2 , − 7 4 2 } }

We can graphically determine the intervals of increase and decrease by noting that if f′(x)>0 (f′(x)<0), a|f′(x)|/f′(x)=a (a|f′(x)|/f′(x)=−a). Similarly, the intervals for which the graph is concave up and concave down can be determined by noting that if f″(x)>0 (f″(x)<0), a|f″(x)|/f″(x)=a (a|f″(x)|/f″(x)=−a). We use Plot to graph |f′(x)|/f′(x) and 2|f″(x)|/f″(x) (different values are used so we can differentiate between the two plots) in Fig. 3.20.

Figure 3.20 Graphs of |f^′(x)|/f^′(x) and 2|f^″(x)|/f^″(x) (Cornell University colors).

Plot [ { Abs [ d1 ] / d1 , 2 Abs [ d2 ] / d2 } , { x , − 2 , 2 } , PlotRange → { − 3 , 3 } ,

PlotStyle→{{CMYKColor[0,1.,.79,.2]},{CMYKColor[.72,.66,.65,72]}}]

From the graph, we see that f′(x)>0 for x in (−∞,−1)∪(1,∞), f′(x)<0 for x in (−1,1), f″(x)>0 for x in (−1/2,0)∪(1/2,∞), and f″(x)<0 for x in (−∞,−1/2)∪(0,1/2). Thus, the graph of f(x) is

•  increasing and concave down for x in (−∞,−1),

•  decreasing and concave down for x in (−1,−1/2),

•  decreasing and concave up for x in (−1/2,0),

•  decreasing and concave down for x in (0,12),

•  decreasing and concave up for x in (1/2,1), and

•  increasing and concave up for x in (1,∞).

We also see that f(0)=0 is neither a relative minimum nor maximum. To see all points of interest, our domain must contain −1 and 1 while our range must contain −2 and 2. We choose to graph f(x) for −2⩽x⩽2; we choose the range displayed to be −4⩽y⩽4. (See Fig. 3.21.)

Figure 3.21 f(x) for −2 ⩽ x ⩽ 2 and −4 ⩽ y ⩽ 4.

Plot [ f [ x ] , { x , − 2 , 2 } , PlotRange → { − 4 , 4 } ,

PlotStyle - > CMYKColor [ 0 , 1 . , . 79 , . 2 ] ,

AxesLabel→{x,y}]  □

Remember to be especially careful when working with functions that involve odd roots. If n is odd and x is even, use Surd[x,n] to obtain the real-valued nth root of x.

Example 3.19

Graph f(x)=(x−2)2/3(x+1)1/3.

Solution

We begin by defining f(x) and then computing and simplifying f′(x) and f″(x) with ′ and Simplify.

Clear [ f ]

f [ x _ ] = Surd [ ( x − 2 ) ∧ 2 , 3 ] Surd [ ( x + 1 ) , 3 ] ;

d1 = Simplify [ f ′ [ x ] ]

d2 = Simplify [ f ” [ x ] ]

( − 2 + x ) x ( − 2 + x ) 2 3 2 1 + x 3 2

− 2 ( 1 + x ) ( − 2 + x ) 2 3 2 1 + x 3 2

By inspection, we see that the critical numbers are x=0, 2, and −1. We cannot use Theorem 3.5 to classify f(2) and f(−1) because f″(x) is undefined if x=2 or −1. On the other hand, f″(0)<0 so f(0)=22/3 is a relative maximum. By hand, we make a sign chart to see that the graph of f(x) is

•  increasing and concave up on (−∞,−1),

•  increasing and concave down on (−1,0),

•  decreasing and concave down on (0,2), and

•  increasing and concave down on (2,∞).

Hence, f(−1)=0 is neither a relative minimum nor maximum while f(2)=0 is a relative minimum by Theorem 3.4. We use Plot to graph f(x) for −2⩽x⩽3 in Fig. 3.22.

Figure 3.22 f(x) for −2 ⩽ x ⩽ 3 (Carnegie Mellon colors).

f [ 0 ]

Plot [ f [ x ] , { x , − 2 , 3 } ,

PlotStyle → CMYKColor [ . 07 , 1 , . 82 , . 26 ] ,

AxesLabel → { x , y } ]

22/3  □

The previous examples illustrate that if x=a is a critical number of f(x) and f′(x) makes a simple change in sign from positive to negative at x=a, then (a,f(a)) is a relative maximum. If f′(x) makes a simple change in sign from negative to positive at x=a, then (a,f(a)) is a relative minimum. Mathematica is especially useful in investigating interesting functions for which this may not be the case.

Example 3.20

Consider

f ( x ) = { x 2 sin 2 ⁡ ( 1 x ) , x ≠ 0 0 , x = 0 ,

x=0 is a critical number because f′(x) does not exist if x=0. The point (0,0) is both a relative and absolute minimum, even though f′(x) does not make a simple change in sign at x=0, as illustrated in Fig. 3.23.

Figure 3.23 f(x)=[xsin⁡(1x)]2 and f^′(x) for −0.1 ⩽ x ⩽ 0.1 (Carnegie Mellon colors).

f [ x _ ] = ( x Sin [ 1 x ] ) 2 ;

f ′ [ x ] // Factor

− 2 Sin [ 1 x ] ( Cos [ 1 x ] − x Sin [ 1 x ] )

p1 = Plot [ f [ x ] , { x , − 0.1 , 0.1 } ,

PlotStyle → CMYKColor [ . 92 , . 02 , 1 , . 12 ] ,

AxesLabel → { x , y } , PlotLabel → “(a)” ] ;

p2 = Plot [ f ′ [ x ] , { x , − 0.1 , 0.1 } ,

PlotStyle → CMYKColor [ 0 , . 32 , 1 , 0 ] ,

AxesLabel → { x , y } , PlotLabel → “(b)” ] ;

Show[GraphicsRow[{p1,p2}]]

Notice that the derivative “oscillates” infinitely many times near x=0, so the first derivative test cannot be used to classify (0,0).

The functions Maximize and Minimize can be used to assist with finding extreme values. For a function of a single variable Maximize[f[x],x] (Minimize[f[x],x]) attempts to find the maximum (minimum) values of f(x); Maximize[{f[x],a<=x<=b},x] (Minimize[{f[x],a<=x<=b},x]) attempts to find the maximum (minimum) values of f(x) on [a,b].

Example 3.21

Consider f(x)=110(−12x+3x2+2x3). After defining f(x), we plot f(x) and f′(x) together in Fig. 3.24.

Figure 3.24 f(x) has one relative maximum and one relative minimum but no absolute extreme values (Purdue University colors).

f [ x _ ] = 1 / 10 ( − 12 x + 3 x ∧ 2 + 2 x ∧ 3 ) ;

Plot [ Tooltip [ { f [ x ] , f ′ [ x ] } ] , { x , − 4 , 4 } , PlotRange → { − 4 , 4 } ,

AspectRatio → Automatic , PlotStyle → { { Black } , { CMYKColor [ . 06 , . 27 , 1 , . 12 ] } } ]

With Maximize, we see that f(x) does not have a maximum on its domain. However, when we restrict the interval to −3⩽x⩽2, Maximize finds the relative maximum at x=−2.

Maximize [ f [ x ] , x ]

Maximize: The maximum is not attained at any point satisfying the given constraints. 

{ ∞ , { x → ∞ } }

Maximize[{f[x],−3⩽x⩽2},x]

{ 2 , { x → − 2 } }

Similarly, with Minimize we see that the f(x) does not have a minimum value on its domain but find the relative minimum when we restrict the interval to −3⩽x⩽2.

Minimize [ f [ x ] , x ]

Minimize: The minimum is not attained at any point satisfying the given constraints.

{ − ∞ , { x → − ∞ } }

Minimize [ { f [ x ] , − 3 ⩽ x ⩽ 2 } , x ]

{ − 7 10 , { x → 1 } }

However, with Solve, we easily find the two zeros of f′(x) that we see in Fig. 3.24

Solve [ f ′ [ x ] = = 0 , x ]

{ { x → − 2 } , { x → 1 } }

When using Maximize or Minimize you should verify your results using another method.

Example 3.22

The function f(x)=x/(x2+1) is continuous on (−∞,∞) and limx→±∞⁡f(x)=0. Thus, f(x) has an absolute minimum and maximum value on its domain. In this case,

Maximize [ x / ( x ∧ 2 + 1 ) , x ]

{ 1 2 , { x → 1 } }

Minimize [ x / ( x ∧ 2 + 1 ) , x ]

{ − 1 2 , { x → − 1 } }

gives us the absolute maximum and minimum values of f(x) and the x-values where they occur. On the other hand, f(x)=x4−x2 is continuous on (−∞,∞) and limx→±∞⁡f(x)=∞. Thus, f(x) has an absolute minimum on its domain. Because the derivative of a fourth degree polynomial is a third degree polynomial, we know that f′(x) has three zeros, two of which probably correspond to relative minimums. Because the graph of f(x) is symmetric with respect to the y-axis, we further suspect that the absolute minimum is obtained twice—at each relative minimum. Maximize and Minimize give us the following results.

A polynomial of degree n has n zeros (counting multiplicity).

Maximize [ x ∧ 4 − x ∧ 2 , x ]

Maximize: The maximum is not attained at any point satisfying the given constraints. 

{∞,{x→−∞}}

Minimize [ x ∧ 4 − x ∧ 2 , x ]

{ − 1 4 , { x → − 1 2 } }

Note that the result returned by Maximize is correct. Similarly, the result returned by Minimize is correct, but a complete answer would indicate that the absolute minimum value occurs at both x=−1/2 and x=1/2.

Example 3.23

The function f(x)=(x+1)2/(x−2) has a vertical asymptote at x=2. From the derivative,

f [ x _ ] = ( x + 1 ) ∧ 2 / ( x − 2 ) ;

d1 = Simplify [ f ′ [ x ] ]

cns = Solve [ f ′ [ x ] = = 0 , x ]

( − 5 + x ) ( 1 + x ) ( − 2 + x ) 2

{ { x → − 1 } , { x → 5 } }

f [ x ] /. cns

{ 0 , 12 }

we find two critical numbers, one of which is a relative maximum and one is a relative minimum. See Fig. 3.25.

Figure 3.25 A function for which a relative minimum has a function value greater than the function value of a relative maximum (University of California–San Diego colors).

Plot [ Tooltip [ { f [ x ] , f ′ [ x ] } ] , { x , − 6 , 10 } ,

PlotStyle → { { CMYKColor [ 1 , . 86 , . 42 , . 42 ] } ,

{CMYKColor[{.06,.35,.99,.18}]}}]

On the other hand, Maximize and Minimize return confusing results because the function is undefined if x=2. The function has relative extreme values but not absolute extreme values.

Maximize [ f [ x ] , x ]

Maximize: The maximum is not attained at any point satisfying the given constraints. 

{ ∞ , { x → 2 } }

Minimize [ f [ x ] , x ]

Minimize: The minimum is not attained at any point satisfying the given constraints.

{ − ∞ , { x → 2 } }

From the graph in Fig. 3.25, we see that limx→2+⁡f(x)=+∞ while limx→2−⁡f(x)=−∞.

For periodic functions, such as sine and cosine, Maximize and Minimize generally don't indicate all extreme values.

Maximize [ Sin [ x ] , x ]

{ 1 , { x → π 2 } }

Minimize [ Cos [ x ] , x ]

{ − 1 , { x → π } }

3.2.6 Applied Max/Min Problems

Mathematica can be used to assist in solving maximization/minimization problems encountered in a differential calculus course.

Example 3.24

A woman is located on one side of a body of water 4 miles wide. Her position is directly across from a point on the other side of the body of water 16 miles from her house, as shown in the following figure.

If she can move across land at a rate of 10 miles per hour and move over water at a rate of 6 miles per hour, find the least amount of time for her to reach her house.

Solution

From the figure, we see that the woman will travel from A to B by land and then from B to D by water. We wish to find the least time for her to complete the trip.

Let x denote the distance BC, where 0⩽x⩽16. Then, the distance AB is given by 16−x and, by the Pythagorean theorem, the distance BD is given by x2+42. Because rate×time=distance, time=distance/rate. Thus, the time to travel from A to B is 110(16−x), the time to travel from B to D is 16x2+16, and the total time to complete the trip, as a function of x, is

t i m e ( x ) = 1 10 ( 16 − x ) + 1 6 x 2 + 16 , 0 ⩽ x ⩽ 16 .

We must minimize the function time. First, we define time and then verify that time has a minimum by graphing time on the interval [0,16] in Fig. 3.26.

Figure 3.26 Plot of time(x)=110(16−x)+16x2+16, 0 ⩽ x ⩽ 16 (The University of Texas at Austin colors).

Clear [ time ]

time [ x _ ] = 16 − x 10 + 1 6 x 2 + 16 ;

Plot [ time [ x ] , { x , 0 , 16 } , PlotRange → { { 0 , 16 } , { 2 , 3 } } ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ 0 , . 65 , 1 , . 09 ] } ,

AxesLabel → { x , t } ]

Next, we compute the derivative of time and find the values of x for which the derivative is 0 with Solve. The resulting output is named critnums using ReplaceAll (\.).

Together [ time ′ [ x ] ]

5 x − 3 16 + x 2 30 16 + x 2

critnums = Solve [ time ′ [ x ] == 0 ]

{ { x → 3 } }

At this point, we can calculate the minimum time by calculating time[3].

time [ 3 ]

32 15

Alternatively, we demonstrate how to find the value of time[x] for the value(s) listed in critnums.

time [ x ] /. x → 3

32 15

Regardless, we see that the minimum time to complete the trip is 32/15 hours. □

One of the more interesting applied max/min problems is the beam problem. We present two solutions.

Example 3.25

The Beam Problem

Find the exact length of the longest beam that can be carried around a corner from a hallway 2 feet wide to a hallway that is 3 feet wide. (See Fig. 3.27.)

Figure 3.27 The length of the beam is found using similar triangles.

Solution

We assume that the beam has negligible thickness. Our first approach is algebraic. Using Fig. 3.27, which is generated with

f [ x _ ] = x + 2 ;

p1 = Plot [ f [ x ] , { x , 0 , 4 } , PlotStyle → { Thickness [ . 01 ] , CMYKColor [ 1 , . 31 , . 08 , . 42 ] } ,

PlotRange - > { 0 , 6 } ] ;

p2 =

Graphics [ { CMYKColor [ . 65 , . 43 , . 26 , . 78 ] ,

Line [ { { 1 , 0 } , { 1 , f [ 1 ] } , { 4 , f [ 1 ] } , { 4 , f [ 4 ] } , { 4 , f [ 4 ] } , { 0 , f [ 4 ] } ,

{0,0},{1,0}}]}];

p3 = Graphics [ { Text [ “2” , { . 5 , . 2 } ] , { Text [ “3” , { 3.8 , 4.5 } ] } } ] ;

p4 = Graphics [ { CMYKColor [ . 62 , . 19 , . 45 , . 5 ] , Dashing [ { 0.01 , 0.01 } ] ,

Line [ { { 0 , f [ 0 ] } , { 1 , f [ 0 ] } } ] } ] ;

p5 = Graphics [ { Text [ “ θ ” , { . 5 , 2.25 } ] , Text [ “ θ ” , { 1.5 , 3.25 } ] } ] ;

p6 = Graphics [ { Text [ x , { . 9 , 2.35 } ] , Text [ y , { 2.5 , 3.25 } ] } ] ;

Show [ p1 , p2 , p3 , p4 , p5 , p6 , Axes - > None ]

and the Pythagorean theorem, the total length of the beam is

L = 2 2 + x 2 + y 2 + 3 2 .

By similar triangles,

y 3 = 2 x so y = 6 x

and the length of the beam, L, becomes

L ( x ) = 4 + x 2 + 9 + 36 x 2 , 0 < x < ∞ .

Observe that the length of the longest beam is obtained by minimizing L. (Why?)

We ignore negative and imaginary values because length must be nonnegative real number.

Clear [ l ] ;

l [ x _ ] = Sqrt [ 2 ∧ 2 + x ∧ 2 ] + Sqrt [ y ∧ 2 + 3 ∧ 2 ] /. y - > 6 / x

9 + 36 x 2 + 4 + x 2

We use two different methods to solve L′(x)=0. Differentiating

l ′ [ x ]

− 36 9 + 36 x 2 x 3 + x 4 + x 2

36 ∧ 2

1296

Solve [ − 12 4 + x 2 + x 4 4 + x 2 x 2 == 0 , x ]

{ { x → − 2 i } , { x → 2 i } , { x → − 2 2 / 3 3 1 / 3 } , { x → 2 2 / 3 3 1 / 3 } }

p1 = x ∧ 8 ( 9 + 36 / x ∧ 2 ) − 1296 ( 4 + x ∧ 2 ) // Expand // Factor

9(4+x2)(−12+x3)(12+x3)

and solving L′(x)=0 gives us

Solve [ p1 == 0 , x ]

{ { x → − 2 i } , { x → 2 i } , { x → − ( − 3 ) 1 / 3 2 2 / 3 } , { x → ( − 3 ) 1 / 3 2 2 / 3 } , { x → − ( − 2 ) 2 / 3 3 1 / 3 } ,

{ x → ( − 2 ) 2 / 3 3 1 / 3 } , { x → − 2 2 / 3 3 1 / 3 } , { x → 2 2 / 3 3 1 / 3 } }

N [ 2 2 / 3 3 1 / 3 ]

2.28943

l [ 2 2 / 3 3 1 / 3 ]

9 + 3 2 2 / 3 3 1 / 3 + 4 + 2 2 1 / 3 3 2 / 3

N [ % ]

7.02348

2 ( 3 / 2 ) ∧ ( 1 / 3 ) Sqrt [ 1 + ( 4 / 9 ) ∧ ( 1 / 3 ) ] + 3 Sqrt [ 1 + ( 4 / 9 ) ∧ ( 1 / 3 ) ] // N

7.02348

It follows that the length of the beam is L(22/331/3)=9+3⋅22/3⋅31/3+4+2⋅21/3⋅32/3=13+9⋅22/3⋅31/3+6⋅21/3⋅32/3≈7.02. See Fig. 3.28.

Figure 3.28 Graph of L(x).

Plot [ l [ x ] , { x , 0 , 20 } , PlotRange - > { 0 , 20 } , AspectRatio - > Automatic ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ . 62 , . 19 , . 45 , . 5 ] } , AxesLabel - > { x , y } ]

Our second approach uses right triangle trigonometry. In terms of θ, the length of the beam is given by

L ( θ ) = 2 csc ⁡ θ + 3 sec ⁡ θ , 0 < θ < π / 2 .

Differentiating gives us

L ′ ( θ ) = − 2 csc ⁡ θ cot ⁡ θ + 3 sec ⁡ θ tan ⁡ θ .

To avoid typing the θ symbol, we define L as a function of t.

Clear [ l ]

l[t_]=2Csc[t]+3Sec[t]

2 Csc [ t ] + 3 Sec [ t ]

We now solve L′(θ)=0. First multiply through by sin⁡θ and then by tan⁡θ.

3 sec ⁡ θ tan ⁡ θ = 2 csc ⁡ θ cot ⁡ θ tan 2 ⁡ θ = 2 3 cot ⁡ θ tan 3 ⁡ θ = 2 3 tan ⁡ θ = 2 3 3 .

In this case, observe that we cannot compute θ exactly. However, we do not need to do so. Let 0<θ<π/2 be the unique solution of tan⁡θ=2/33. See Fig. 3.29. Using the identity tan2⁡θ+1=sec2⁡θ, we find that sec⁡θ=1+4/93. Similarly, because cot⁡θ=3/23 and cot2⁡θ+1=csc2⁡θ, csc⁡θ=3/231+4/93. Hence, the length of the beam is

L ( θ ) = 2 3 2 3 1 + 4 9 3 + 3 1 + 4 9 3 ≈ 7.02 .

Figure 3.29 Graph of L(θ) and L^′(θ).

Plot [ Tooltip [ { l [ t ] , l ′ [ t ] } ] , { t , 0 , Pi / 2 } , PlotRange - > { − 20 , 20 } ,

PlotStyle - > { { Thickness [ . 01 ] , CMYKColor [ . 52 , . 59 , . 45 , . 9 ] } ,

{ Thickness [ . 01 ] , CMYKColor [ . 03 , . 04 , . 14 , . 08 ] } } ,

When you use Tooltip, scrolling the cursor over the plot will identify the plot for you.

AxesLabel→{t,y}]  □

In the next two examples, the constants do not have specific numerical values.

Example 3.26

Find the volume of the right circular cone of maximum volume that can be inscribed in a sphere of radius R.

Solution

Try to avoid three-dimensional figures unless they are absolutely necessary. For this problem, a cross-section of the situation is sufficient. See Fig. 3.30, which is created with

Figure 3.30 Cross-section of a right circular cone inscribed in a sphere (UCLA colors).

p1 = ParametricPlot [ { Cos [ t ] , Sin [ t ] } , { t , 0 , 2 Pi } ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 0 , . 05 , 1 , 0 ] } } ] ;

p2 = Graphics [ { CMYKColor [ . 75 , . 35 , 0 , . 07 ] , Thickness [ . 005 ] ,

Line [ { { 0 , 1 } , { Cos [ 4 Pi / 3 ] , Sin [ 4 Pi / 3 ] } , { Cos [ 5 Pi / 3 ] , Sin [ 5 Pi / 3 ] } , { 0 , 1 } } ] ,

PointSize [ . 02 ] , Point [ { 0 , 0 } ] , Line [ { { Cos [ 4 Pi / 3 ] , Sin [ 4 Pi / 3 ] } , { 0 , 0 } , { 0 , 1 } } ] ,

Line [ { { 0 , 0 } , { 0 , Sin [ 4 Pi / 3 ] } } ] } ] ;

p3 = Graphics [ { Text [ R , { − . 256 , − . 28 } ] , Text [ R , { − . 04 , . 5 } ] , Text [ y , { − . 04 , − . 5 } ] ,

Text [ x , { − . 2 , − . 8 } ] } ] ;

Show [ p1 , p2 , p3 , AspectRatio - > Automatic , Ticks - > None , Axes - > None ]

The volume, V, of a right circular cone with radius r and height h is V=13πr2h. Using the notation in Fig. 3.30, the volume is given by

V = 1 3 π x 2 ( R + y ) .

(3.4)

However, by the Pythagorean theorem, x2+y2=R2 so x2=R2−y2 and equation (3.4) becomes

V = 1 3 π ( R 2 − y 2 ) ( R + y ) = 1 3 π ( R 3 + R 2 y − R y 2 − y 3 ) ,

(3.5)

s1 = Expand [ ( r ∧ 2 − y ∧ 2 ) ( r + y ) ]

r3+r2y−ry2−y3

where 0⩽y⩽R. V(y) is continuous on [0,R] so it will have a minimum and maximum value on this interval. Moreover, the minimum and maximum values either occur at the endpoints of the interval or at the critical numbers on the interior of the interval. Differentiating equation (3.5) with respect to y gives us

Remember that R>0 is a constant.

d V d y = 1 3 π ( R 2 − 2 R y − 3 y 2 ) = 1 3 π ( R − 3 y ) ( R + y )

s2 = D [ s1 , y ]

r 2 − 2 r y − 3 y 2

and we see that dV/dy=0 if y=13R or y=−R.

Factor [ s2 ]

( r − 3 y ) ( r + y )

Solve [ s2 = = 0 , y ]

{ { y → − r } , { y → r 3 } }

We ignore y=−R because −R is not in the interval [0,R]. Note that V(0)=V(R)=0. The maximum volume of the cone is

V ( 1 3 R ) = 1 3 π ⋅ 32 27 R 3 = 32 81 π R 2 ≈ 1.24 R 3 .

s3 = s1 /. y - > r / 3 // Together

32r327

s3 ⁎ 1 / 3 Pi

32 π r 3 81

N [ % ]

1.24112r3  □

Example 3.27

The Stayed-Wire Problem

Two poles D feet apart with heights L1 feet and L2 feet are to be stayed by a wire as shown in Fig. 3.31. Find the minimum amount of wire required to stay the poles, as illustrated in Fig. 3.31, which is generated with

Figure 3.31 When the wire is stayed to minimize the length, the result is two similar triangles (Texas A & M colors).

p1 = Graphics [ { CMYKColor [ . 15 , 1 , . 39 , . 69 ] , Thickness [ . 0075 ] ,

Line [ { { 0 , 0 } , { 0 , 4 } , { 3.5 , 0 } , { 9 , 5.5 } , { 9 , 0 } , { 0 , 0 } } ] } ] ;

p2 = Graphics [ { Text [ " L 1 " , { . 2 , 2 } ] , Text [ " L 2 " , { 8.8 , 2.75 } ] , Text [ x , { 1.75 , . 2 } ] ,

Text [ x , { 1.75 , . 2 } ] , Text [ " L 1 2 + x 2 " , { 1.75 , 2.75 } ] ,

Text [ " ( D − x ) 2 + L 2 2 " , { 5.5 , 2.75 } ] , Text [ “D-x” , { 6.5 , . 2 } ] } ] ;

Show [ p1 , p2 ]

Solution

Using the notation in Fig. 3.31, the length of the wire, L, is

L ( x ) = L 1 2 + x 2 + L 2 2 + ( D − x ) 2 , 0 ⩽ x ⩽ D .

(3.6)

In the special case that L1=L2, the length of the wire to stay the beams is minimized when the wire is placed halfway between the two beams, at a distance D/2 from each beam. Thus, we assume that the lengths of the beams are different; we assume that L1<L2, as illustrated in Fig. 3.31. We compute L′(x) and then solve L′(x)=0. We use PowerExpand because PowerExpand[expr] expands out all products and powers assuming the variables are real and positive. That is, with PowerExpand we obtain that x2=x rather than x2=|x|.

Clear [ l ]

l [ x _ ] = Sqrt [ x ∧ 2 + l1 ∧ 2 ] + Sqrt [ ( d − x ) ∧ 2 + l1 ∧ 2 ]

l1 2 + ( d − x ) 2 + l1 2 + x 2

l ′ [ x ] // Together

− d l1 2 + x 2 + x l1 2 + x 2 + x d 2 + l1 2 − 2 d x + x 2 l1 2 + x 2 d 2 + l1 2 − 2 d x + x 2

Solve [ l ′ [ x ] == 0 , x ]

{ { x → d 2 } }

The result indicates that x=L1D/(L1+L2) minimizes L(x). (Note that we ignore the other value because L1−L2<0.) Moreover, the triangles formed by minimizing L are similar triangles.

Clear [ l ]

l [ x _ ] = Sqrt [ x ∧ 2 + l1 ∧ 2 ] + Sqrt [ ( d − x ) ∧ 2 + l1 ∧ 2 ]

l1 2 + ( d − x ) 2 + l1 2 + x 2

l ′ [ x ] // Together

− d l1 2 + x 2 + x l1 2 + x 2 + x d 2 + l1 2 − 2 d x + x 2 l1 2 + x 2 d 2 + l1 2 − 2 d x + x 2

l [ 0 ] // PowerExpand

l1 + d 2 + l1 2

l [ d ] // PowerExpand

l1 + d 2 + l1 2

Solve [ l ′ [ x ] == 0 , x ]

{ { x → d 2 } }

l [ d / 2 ] // Together // PowerExpand

d 2 + 4 l1 2

d 2 + 4 l1 2

d2+4l12

l1 / ( d l1 l1 + l2 ) // Simplify

l1 + l2 d

l2 / ( d − d l1 l1 + l2 ) // Simplify

l1+l2d  □

3.2.7 Antidifferentiation

3.2.7.1 Antiderivatives

F(x) is an antiderivative of f(x) if F′(x)=f(x). The symbols

∫ f ( x ) d x

mean “find all antiderivatives of f(x).” Because all antiderivatives of a given function differ by a constant, we usually find an antiderivative, F(x), of f(x) and then write

∫ f ( x ) d x = F ( x ) + C ,

where C represents an arbitrary constant. The command

Integrate[f[x],x]

attempts to find an antiderivative, F(x), of f(x). Instead of using Integrate, you might prefer to use the button on the Basic Math Input or Basic Math Assistant palettes to help you evaluate antiderivatives. Mathematica does not include the “+C” that we include when writing ∫f(x)dx=F(x)+C. In the same way as D can differentiate many functions, Integrate can antidifferentiate many functions. However, antidifferentiation is a fundamentally difficult procedure so it is not difficult to find functions f(x) for which the command Integrate[f[x],x] returns unevaluated.

Example 3.28

Evaluate each of the following antiderivatives: (a) ∫1x2e1/xdx, (b) ∫x2cos⁡xdx, (c) ∫x21+x2dx, (d) ∫x2−x+2x3−x2+x−1dx, and (e) ∫sin⁡xxdx.

Solution

Entering

Integrate [ 1 / x ∧ 2 Exp [ 1 / x ] , x ]

− e 1 x

shows us that ∫1x2e1/xdx=−e1/x+C. To use the button, first click on the button, fill in the blanks, and press Enter.

Notice that Mathematica does not automatically include the arbitrary constant, C. When computing several antiderivatives, you can use Map to apply Integrate to a list of antiderivatives. However, because Integrate is threadable,

Map[Integrate[#,x]&,list]

returns the same result as Integrate[list,x], which we illustrate to compute (b), (c), and (d).

Integrate [ { x ∧ 2 Cos [ x ] , x ∧ 2 Sqrt [ 1 + x ∧ 2 ] ,

( x ∧ 2 − x + 2 ) / ( x ∧ 3 − x ∧ 2 + x − 1 ) } , x ]

{ 2 x Cos [ x ] + ( − 2 + x 2 ) Sin [ x ] , 1 8 ( 1 + x 2 ( x + 2 x 3 ) − ArcSinh [ x ] ) ,

− ArcTan [ x ] + Log [ 1 − x ] }

For (e), we see that there is not a “closed form” antiderivative of ∫sin⁡xxdx and the result is given in terms of a definite integral, the sine integral function:

S i ( x ) = ∫ 0 x sin ⁡ t t d t .

Integrate [ Sin [ x ] / x , x ]

SinIntegral[x]  □

u-Substitutions

Usually, the first antidifferentiation technique discussed is the method of u-substitution. Suppose that F(x) is an antiderivative of f(x). Given

∫ f ( g ( x ) ) g ′ ( x ) d x ,

we let u=g(x) so that du=g′(x)dx. Then,

∫ f ( g ( x ) ) g ′ ( x ) d x = ∫ f ( u ) d u = F ( u ) + C = F ( g ( x ) ) + C ,

where F(x) is an antiderivative of f(x). After mastering u-substitutions, the integration by parts formula,

∫ u d v = u v − ∫ v d u ,

(3.7)

is introduced.

Example 3.29

Evaluate ∫2x4x−1dx.

Solution

We use Integrate to evaluate the antiderivative.

i1 = Integrate [ 2 ∧ x Sqrt [ 4 ∧ x − 1 ] , x ]

2 x − 1 + 4 x − Log [ 2 x + − 1 + 4 x ] Log [ 4 ]

Proceeding by hand, we let u=2x. Then, du=2xln⁡2dx or, equivalently, 1ln⁡2du=2xdx

D [ 2 ∧ x , x ]

2 x Log [ 2 ]

so ∫2x4x−1dx=1ln⁡2∫u2−1du. We now use Integrate to evaluate 1ln⁡2∫u2−1du

i2 = 1 / Log [ 2 ] Integrate [ Sqrt [ u ∧ 2 − 1 ] , u ]

1 2 u − 1 + u 2 − 1 2 Log [ u + − 1 + u 2 ] Log [ 2 ]

Simplify [ i2 ]

u − 1 + u 2 − Log [ u + − 1 + u 2 ] Log [ 4 ]

and then /. (ReplaceAll)/ to replace u with 2x.

i3 = i2 /. u → 2 ∧ x

2 − 1 + x − 1 + 2 2 x − 1 2 Log [ 2 x + − 1 + 2 2 x ] Log [ 2 ]

Observe that the result we obtained by hand is the same as the result obtained by Integrate directly. Sometimes, the results will look different and have slightly different forms. To verify that they are equivalent, subtract the two, and simplify the result. If the result is a constant, the two antiderivatives are equivalent. If not, they aren't.  □

As we did with derivatives, with DynamicModule, we create a simple dynamic that lets you compute the derivative and antiderivative of basic functions and plot them on a standard viewing window, [−5,5]×[−5,5]. The layout of Fig. 3.32 is primarily determined by Panel, Column, and Grid.

Figure 3.32 Seeing the relationship between the derivative and antiderivative of a function and the original function (Texas A & M colors).

Panel [ DynamicModule [ { f = x ∧ 2 } ,

Column [ { InputField [ Dynamic [ f ] ] , Grid [ { { “First Derivative” ,

Panel [ Dynamic [ D [ f , x ] // Simplify ] ] } ,

{ “Antiderivative” ,

Panel [ Dynamic [ Integrate [ f , x ] // Simplify ] ] } } ] ,

Dynamic [ Plot [ Evaluate [ Tooltip [ { f , D [ f , x ] ,

Integrate [ f , x ] } ] ] , { x , − 5 , 5 } , PlotRange → { − 5 , 5 } ,

AspectRatio → Automatic , AxesLabel → { x , y } ,

PlotStyle → { { CMYKColor [ . 15 , 1 , . 39 , . 69 ] } , { CMYKColor [ 1 , . 48 , . 09 , . 46 ] } ,

{ CMYKColor [ . 46 , . 23 , . 84 , . 68 ] } } ] ] } ] ] , ImageSize → { 300 , 300 } ]

3.3 Integral Calculus

3.3.1 Area

In integral calculus courses, the definite integral is frequently motivated by investigating the area under the graph of a positive continuous function on a closed interval. Let y=f(x) be a nonnegative continuous function on an interval [a,b] and let n be a positive integer. If we divide [a,b] into n subintervals of equal length and let [xk−1,xk] denote the kth subinterval, the length of each subinterval is (b−a)/n and xk=a+kb−an. The area bounded by the graphs of y=f(x), x=a, x=b, and the y-axis can be approximated with the sum

∑ k = 1 n f ( x k ⁎ ) b − a n ,

(3.8)

where xk⁎∈[xk−1,xk]. Typically, we take xk⁎=xk−1=a+(k−1)b−an (the left endpoint of the kth subinterval), xk⁎=xk−1=a+kb−an (the right endpoint of the kth subinterval), or xk⁎=12(xk−1+xk)=a+12(2k−1)b−an (the midpoint of the kth subinterval). For these choices of xk⁎, (3.8) becomes

b − a n ∑ k = 1 n f ( a + ( k − 1 ) b − a n )

(3.9)

b − a n ∑ k = 1 n f ( a + k b − a n ) ,  and 

(3.10)

b − a n ∑ k = 1 n f ( a + 1 2 ( 2 k − 1 ) b − a n ) ,

(3.11)

respectively. If y=f(x) is increasing on [a,b], (3.9) is an under approximation and (3.10) is an upper approximation: (3.9) corresponds to an approximation of the area using n inscribed rectangles; (3.10) corresponds to an approximation of the area using n circumscribed rectangles. If y=f(x) is decreasing on [a,b], (3.10) is an under approximation and (3.9) is an upper approximation: (3.10) corresponds to an approximation of the area using n inscribed rectangles; (3.9) corresponds to an approximation of the area using n circumscribed rectangles.

In the following example, we define the functions leftsum[f[x],a,b,n], middlesum[f[x],a,b,n], and rightsum[f[x],a,b,n] to compute (3.9), (3.11), and (3.10), respectively, and leftbox[f[x],a,b,n], middlebox[f[x], a,b,n], and rightbox[f[x],a, b,n] to generate the corresponding graphs. After you have defined these functions, you can use them with functions y=f(x) that you define.

Remark 3.6

To define a function of a single variable, f(x)=expressioninx, enter f[x_]=expression in x. To generate a basic plot of y=f(x) for a⩽x⩽b, enter Plot[f[x],{x,a,b}].

Example 3.30

Let f(x)=9−4x2. Approximate the area bounded by the graph of y=f(x), x=0, x=3/2, and the y-axis using (a) 100 inscribed and (b) 100 circumscribed rectangles. (c) What is the exact value of the area?

Solution

We begin by defining and graphing y=f(x) in Fig. 3.33.

Figure 3.33 f(x) for 0 ⩽ x ⩽ 3/2 (Princeton University colors).

f [ x _ ] = 9 − 4 x ∧ 2 ;

Plot [ f [ x ] , { x , 0 , 3 / 2 } ,

AxesLabel → { x , y } ,

PlotStyle → { { CMYKColor [ 0 , . 61 , . 97 , 0 ] , Thickness [ . 01 ] } } ]

The first derivative, f′(x)=−8x is negative on the interval so f(x) is decreasing on [0,3/2]. Thus, an approximation of the area using 100 inscribed rectangles is given by (3.10) while an approximation of the area using 100 circumscribed rectangles is given by (3.9). After defining leftsum, rightsum, and middlesum, these values are computed using leftsum and rightsum. The use of middlesum is illustrated as well. Approximations of the sums are obtained with N.

N[number] returns a numerical approximation of number.

leftsum [ f _ , a _ , b _ , n _ ] := Module [ { } ,

( b − a ) / n Sum [ f /. x - > a + ( k − 1 ) ( b − a ) / n , { k , 1 , n } ] ] ;

rightsum [ f _ , a _ , b _ , n _ ] := Module [ { } ,

( b − a ) / n Sum [ f /. x - > a + k ( b − a ) / n , { k , 1 , n } ] ] ;

middlesum [ f _ , a _ , b _ , n _ ] := Module [ { } ,

( b − a ) / n Sum [ f /. x - > a + 1 / 2 ( 2 k − 1 ) ( b − a ) / n , { k , 1 , n } ] ] ;

l100 = leftsum [ f [ x ] , 0 , 3 / 2 , 100 ]

N [ % ]

r100 = rightsum [ f [ x ] , 0 , 3 / 2 , 100 ]

N [ % ]

m100 = middlesum [ f [ x ] , 0 , 3 / 2 , 100 ]

N [ % ]

362691 40000

9.06728

357291 40000

8.93228

720009 80000

9.00011

Observe that these three values appear to be close to 9. In fact, 9 is the exact value of the area of the region bounded by y=f(x), x=0, x=3/2, and the y-axis. To help us see why this is true, we define leftbox, middlebox, and rightbox, and then use these functions to visualize the situation using n=4, 16, and 32 rectangles in Fig. 3.34.

It is not important that you understand the syntax of these three functions at this time. Once you have entered the code, you can use them to visualize the process for your own functions, y=f(x).

Figure 3.34 f(x) with 4, 16, and 32 rectangles.

leftbox [ f _ , a _ , b _ , n _ , opts _ _ _ ] := Module [ { z , p1 , recs , ls } ,

z [ k _ ] = a + ( b − a ) k / n ;

p1 = Plot [ f , { x , a , b } , PlotRange → All ,

PlotStyle - > { { Thickness [ . 02 ] , CMYKColor [ 0 , . 61 , . 97 , 0 ] } } ] ;

recs = Table [ Rectangle [ { z [ k − 1 ] , 0 } , { z [ k ] , f /. x - > z [ k − 1 ] } ] , { k , 1 , n } ] ;

ls = Table [ Line [ { { z [ k − 1 ] , 0 } , { z [ k − 1 ] , f /. x - > z [ k − 1 ] } , { z [ k ] , f /. x - > z [ k − 1 ] } ,

{ z [ k ] , 0 } } ] , { k , 1 , n } ] ;

Show [ Graphics [ { Black , recs } ] , Graphics [ ls ] , p1 , opts , Axes - > Automatic ,

AspectRatio → 1 ] ]

rightbox [ f _ , a _ , b _ , n _ , opts _ _ _ ] := Module [ { z , p1 , recs , ls } ,

z [ k _ ] = a + ( b − a ) k / n ;

p1 = Plot [ f , { x , a , b } , PlotRange → All ,

PlotStyle - > { { Thickness [ . 02 ] , CMYKColor [ 0 , . 61 , . 97 , 0 ] } } ] ;

recs = Table [ Rectangle [ { z [ k − 1 ] , 0 } , { z [ k ] , f /. x - > z [ k ] } ] , { k , 1 , n } ] ;

ls = Table [ Line [ { { z [ k − 1 ] , 0 } , { z [ k − 1 ] , f /. x - > z [ k ] } , { z [ k ] , f /. x - > z [ k ] } ,

{ z [ k ] , 0 } } ] , { k , 1 , n } ] ;

Show [ Graphics [ { Black , recs } ] , Graphics [ ls ] , p1 , opts ,

Axes - > Automatic , AspectRatio → 1 ] ]

middlebox [ f _ , a _ , b _ , n _ , opts _ _ _ ] := Module [ { z , p1 , recs , ls } ,

z [ k _ ] = a + ( b − a ) k / n ;

p1 = Plot [ f , { x , a , b } , PlotRange → All ,

PlotStyle - > { { Thickness [ . 02 ] , CMYKColor [ 0 , . 61 , . 97 , 0 ] } } ] ;

recs = Table [ Rectangle [ { z [ k − 1 ] , 0 } , { z [ k ] , f /. x - > 1 / 2 ( z [ k − 1 ] + z [ k ] ) } ] ,

{ k , 1 , n } ] ;

ls = Table [ Line [ { { z [ k − 1 ] , 0 } , { z [ k − 1 ] , f /. x - > 1 / 2 ( z [ k − 1 ] + z [ k ] ) } ,

{ z [ k ] , f /. x - > 1 / 2 ( z [ k − 1 ] + z [ k ] ) } , { z [ k ] , 0 } } ] , { k , 1 , n } ] ;

Show [ Graphics [ { Black , recs } ] , Graphics [ ls ] , p1 , opts ,

Axes - > Automatic , AspectRatio → 1 ] ]

somegraphs = { { leftbox [ f [ x ] , 0 , 3 2 , 4 , DisplayFunction → Identity ] ,

middlebox [ f [ x ] , 0 , 3 2 , 4 , DisplayFunction → Identity ] ,

rightbox [ f [ x ] , 0 , 3 2 , 4 , DisplayFunction → Identity ] } ,

{ leftbox [ f [ x ] , 0 , 3 2 , 16 , DisplayFunction → Identity ] ,

middlebox [ f [ x ] , 0 , 3 2 , 16 , DisplayFunction → Identity ] ,

rightbox [ f [ x ] , 0 , 3 2 , 16 , DisplayFunction → Identity ] } ,

{ leftbox [ f [ x ] , 0 , 3 2 , 32 , DisplayFunction → Identity ] ,

middlebox [ f [ x ] , 0 , 3 2 , 32 , DisplayFunction → Identity ] ,

rightbox [ f [ x ] , 0 , 3 2 , 32 , DisplayFunction → Identity ] } } ;

Show [ GraphicsGrid [ somegraphs ] ]

somegraphs = { { leftbox [ f [ x ] , 0 , 3 2 , 4 ] ,

middlebox [ f [ x ] , 0 , 3 2 , 4 ] ,

rightbox [ f [ x ] , 0 , 3 2 , 4 ] } ,

{ leftbox [ f [ x ] , 0 , 3 2 , 16 ] ,

middlebox [ f [ x ] , 0 , 3 2 , 16 ] ,

rightbox [ f [ x ] , 0 , 3 2 , 16 ] } ,

{ leftbox [ f [ x ] , 0 , 3 2 , 32 ] ,

middlebox [ f [ x ] , 0 , 3 2 , 32 ] ,

rightbox [ f [ x ] , 0 , 3 2 , 32 ] } } ;

Show [ GraphicsGrid [ somegraphs ] ]

Notice that as n increases, the under approximations increase to the value of the area while the upper approximations decrease to the value of the area. In the limit as n→∞, if the two limits are equal, we can conclude that the area is the value of the limits.

These graphs help convince us that the limit of the sum as n→∞ of the areas of the inscribed and circumscribed rectangles is the same. We compute the exact value of (3.9) with leftsum, evaluate and simplify the sum with Simplify, and compute the limit as n→∞ with Limit. We see that the limit is 9.

ls = leftsum [ f [ x ] , 0 , 3 / 2 , n ]

ls2 = Simplify [ ls ]

Limit [ ls2 , n - > Infinity ]

9 ( − 1 + 3 n + 4 n 2 ) 4 n 2

9 ( − 1 + 3 n + 4 n 2 ) 4 n 2

9

Similar calculations are carried out for (3.10) and again we see that the limit is 9. We conclude that the exact value of the area is 9.

rs = rightsum [ f [ x ] , 0 , 3 / 2 , n ]

rs2 = Simplify [ rs ]

Limit [ rs2 , n - > Infinity ]

9 ( − 1 − 3 n + 4 n 2 ) 4 n 2

9 ( − 1 − 3 n + 4 n 2 ) 4 n 2

9

For illustrative purposes, we confirm this result with middlesum.

ms = middlesum [ f [ x ] , 0 , 3 / 2 , n ]

ms2 = Simplify [ ms ]

Limit [ ms2 , n - > Infinity ]

9 ( 1 + 8 n 2 ) 8 n 2

9 + 9 8 n 2

9

As illustrated earlier, with Manipulate, you can experiment with different functions and different n values. First, we define a set of “typical” functions.

quad [ x _ ] = 100 − x ∧ 2 ;

cubic [ x _ ] = 4 / 9 x ∧ 3 − 49 / 9 x ∧ 2 + 100 ;

rational [ x _ ] = 100 / ( x ∧ 2 + 1 ) ;

root [ x _ ] = Sqrt [ 10 − x ] ;

sin⁡[x_]=75Sin[Pix/5];

Next, we use Manipulate to create an object that allows us to experiment with how “typical” functions react to changes in n using left, middle, and right-hand endpoint approximations for computations of Riemann sums. In the resulting Manipulate object, n=4 rectangles is the default; you can choose n-values from 0 to 100. The value of the corresponding Riemann sum is shown below the graphic. See Fig. 3.35.

Figure 3.35 With Manipulate, we can investigate Riemann sum approximations and their graphical representations for various functions.

How does the Manipulate object change if you remove Transpose from the command?

Manipulate [ Show [ GraphicsGrid [ { { leftbox [ f [ x ] , 0 , 10 , n ] ,

Graphics [ { Inset [ leftsum [ f [ x ] , 0 , 10 , n ] // N , { 0 , 0 } ] } ] } ,

{ middlebox [ f [ x ] , 0 , 10 , n ] ,

Graphics [ { Inset [ middlesum [ f [ x ] , 0 , 10 , n ] // N , { 0 , 0 } ] } ] } ,

{ rightbox [ f [ x ] , 0 , 10 , n ] ,

Graphics [ { Inset [ rightsum [ f [ x ] , 0 , 10 , n ] // N , { 0 , 0 } ] } ] } } //

Transpose ] ] , { { f , quad } , { quad , cubic , rational , root , sin ⁡ } } ,

{{n,4},0,100,1}]  □

3.3.2 The Definite Integral

In integral calculus courses, we formally learn that the definite integral of the function y=f(x) from x=a to x=b is

∫ a b f ( x ) d x = lim | P | → 0 ⁡ ∑ k = 1 n f ( x k ⁎ ) Δ x k ,

(3.12)

provided that the limit exists. In equation (3.12), P={a=x0<x1<x2<…<xn=b} is a partition of [a,b], |P| is the norm of P,

| P | = max ⁡ { x k − x k − 1 | k = 1 , 2 , … , n } ,

Δxk=xk−xk−1, and xk⁎∈[xk−1,xk].

The Fundamental Theorem of Calculus provides the fundamental relationship between differentiation and integration.

Theorem 3.6

The Fundamental Theorem of Calculus

Assume that y=f(x) is continuous on [a,b].

1.  If F(x)=∫axf(t)dt, then F(x) is an antiderivative of f(x): F′(x)=f(x) or, equivalently,

d d x ( ∫ a x f ( t ) d t ) = f ( x ) .

2.  If G(x) is any antiderivative of f(x), then ∫abf(x)dx=G(b)−G(a).

By the Fundamental Theorem of Calculus and the Chain Rule, it follows that

d d x ( ∫ g ( x ) h ( x ) f ( t ) d t ) = f ( h ( x ) ) ⋅ h ′ ( x ) − f ( g ( x ) ) ⋅ g ′ ( x ) .

Mathematica's Integrate command can compute many definite integrals. The command

Integrate[f[x],{x,a,b}]

attempts to compute ∫abf(x)dx while

Integrate[f[x],x]

attempts to find an antiderivative of f(x). Remember that Mathematica does not include the “+C” when computing antiderivatives. Because integration is a fundamentally difficult procedure, it is easy to create integrals for which the exact value cannot be found explicitly. In those cases, use N to obtain an approximation of its value or obtain a numerical approximation of the integral directly with

NIntegrate[f[x],{x,a,b}].

In the same way as you use the button to compute antiderivatives, you can use the button to compute definite integrals. If the result returned is unevaluated, use N to obtain a numerical approximation of the value of the integral or use NIntegrate.

Example 3.31

Evaluate (a) ∫14(x2+1)/xdx; (b) ∫0π/2xcos⁡x2dx; (c) ∫0πe2xsin2⁡2xdx; (d) ∫012πe−x2dx; and (e) ∫−10u3du.

Solution

We evaluate (a)–(c) directly with Integrate.

Integrate [ ( x ∧ 2 + 1 ) / Sqrt [ x ] , { x , 1 , 4 } ]

72 5

Integrate [ x Cos [ x ∧ 2 ] , { x , 0 , Sqrt [ Pi / 2 ] } ]

1 2

Integrate [ Exp [ 2 x ] Sin [ 2 x ] ∧ 2 , { x , 0 , Pi } ]

15(−1+e2π)

For (d), the result returned is in terms of the error function, Erf[x], which is defined by the integral

E r f [ x ] = 2 π ∫ 0 x e − t 2 d t .

Integrate [ 2 / Sqrt [ Pi ] Exp [ − x ∧ 2 ] , { x , 0 , 1 } ]

Erf [ 1 ]

We use N to obtain an approximation of the value of the definite integral.

Integrate [ 2 / Sqrt [ Pi ] Exp [ − x ∧ 2 ] , { x , 0 , 1 } ] // N

0.842701

(e) Recall that Mathematica does not return a real number when we compute odd roots of negative numbers so the following result would be surprising to many students in an introductory calculus course because it is complex.

Integrate [ u ∧ ( 1 / 3 ) , { u , − 1 , 0 } ]

3 4 ( − 1 ) 1 / 3

Use Surd[u,n] to return the real nth root of u, un, if n is odd.

Integrate [ Surd [ u , 3 ] , { u , − 1 , 0 } ]

−34  □

Improper integrals are computed using Integrate in the same way as with definite integrals.

Example 3.32

Evaluate (a) ∫01ln⁡xxdx; (b) ∫−∞∞2πe−x2dx; (c) ∫1∞1xx2−1dx; (d) ∫0∞1x2+x4dx; (e) ∫241(x−3)23dx; and (f) ∫−∞∞1x2+x−6dx.

Solution

(a) This is an improper integral because the integrand is discontinuous on the interval [0,1] but we see that the improper integral converges to −4.

Integrate [ Log [ x ] / Sqrt [ x ] , { x , 0 , 1 } ]

−4

(b) This is an improper integral because the interval of integration is infinite but we see that the improper integral converges to 2.

Integrate [ 2 / Sqrt [ Pi ] Exp [ − x ∧ 2 ] , { x , − Infinity , Infinity } ]

2

(c) This is an improper integral because the integrand is discontinuous on the interval of integration and because the interval of integration is infinite but we see that the improper integral converges to π/2.

Integrate [ 1 / ( x Sqrt [ x ∧ 2 − 1 ] ) , { x , 1 , Infinity } ]

π 2

(d) As with (c), this is an improper integral because the integrand is discontinuous on the interval of integration and because the interval of integration is infinite but we see that the improper integral diverges to ∞.

(e) Recall that Mathematica does not return a real number when we compute odd roots of negative numbers so the following result would be surprising to many students in an introductory calculus course because it contains imaginary numbers.

Integrate [ 1 / ( x − 3 ) ∧ ( 2 / 3 ) , { x , 2 , 4 } ]

3 − 3 ( − 1 ) 1 / 3

Therefore, we use Surd to obtain the real-valued third root of x−3.

Integrate [ Surd [ 1 / ( x − 3 ) ∧ 2 , 3 ] , { x , 2 , 4 } ]

6

(f) In this case, Mathematica warns us that the improper integral diverges.

To help us understand why the improper integral diverges, we note that 1x2+x−6=15(1x−2−1x+3) and

∫ 1 x 2 + x − 6 d x = ∫ 1 5 ( 1 x − 2 − 1 x + 3 ) d x = 1 5 ln ⁡ ( x − 2 x + 3 ) + C .

Integrate [ 1 / ( x ∧ 2 + x − 6 ) , x ]

1 5 Log [ 2 − x ] − 1 5 Log [ 3 + x ]

Hence the integral is improper because the interval of integration is infinite and because the integrand is discontinuous on the interval of integration so

∫ − ∞ ∞ 1 x 2 + x − 6 d x = ∫ − ∞ − 4 1 x 2 + x − 6 d x + ∫ − 4 − 3 1 x 2 + x − 6 d x + ∫ − 3 0 1 x 2 + x − 6 d x + ∫ 0 2 1 x 2 + x − 6 d x + ∫ 2 3 1 x 2 + x − 6 d x + ∫ 3 ∞ 1 x 2 + x − 6 d x

(3.13)

Evaluating each of these integrals,

we conclude that the improper integral diverges because at least one of the improper integrals in (3.13) diverges. □

In many cases, Mathematica can help illustrate the steps carried out when computing integrals using standard methods of integration like u-substitutions and integration by parts.

Example 3.33

Evaluate (a) ∫ee31xln⁡xdx and (b) ∫0π/4xsin⁡2xdx.

Solution

(a) We let u=ln⁡x. Then, du=1/xdx so

∫ e e 3 1 x ln ⁡ x d x = ∫ 1 3 1 u d u = ∫ 1 3 u − 1 / 2 d u ,

which we evaluate with Integrate.

The new lower limit of integration is 1 because if x=e, u=ln⁡e=1. The new upper limit of integration is 3 because if x=e3, u=ln⁡e3=3.

Integrate [ 1 / Sqrt [ u ] , { u , 1 , 3 } ]

2 ( − 1 + 3 )

To evaluate (b), we let u=x⇒du=dx and dv=sin⁡2xdx⇒v=−12cos⁡2x.

u = x ;

dv = Sin [ 2 x ] ;

du = D [ x , x ]

v = Integrate [ Sin [ 2 x ] , x ]

1

− 1 2 Cos [ 2 x ]

The results mean that

∫ 0 π / 4 x sin ⁡ 2 x d x = − 1 2 x cos ⁡ 2 x ] 0 π / 4 + 1 2 ∫ 0 π / 4 cos ⁡ 2 x d x = 0 + 1 2 ∫ 0 π / 4 cos ⁡ 2 x d x .

The resulting indefinite integral is evaluated with Integrate.

Integrate [ x Sin [ 2 x ] , x ]

− 1 2 x Cos [ 2 x ] + 1 4 Sin [ 2 x ]

Alternatively, we can illustrate the integration by parts calculation, ∫udv=uv−v∫du.

u v − Integrate [ v du , x ]

− 1 2 x Cos [ 2 x ] + 1 4 Sin [ 2 x ]

We use Integrate to evaluate the definite integral.

Integrate [ x Sin [ 2 x ] , { x , 0 , Pi / 4 } ]

14  □

3.3.3 Approximating Definite Integrals

Because integration is a fundamentally difficult procedure to produce exact answers or results, Mathematica is unable to compute a “closed form” of the value of many definite integrals. In these cases, numerical integration can be used to obtain an approximation of the definite integral using N together with Integrate or NIntegrate:

NIntegrate[f[x],{x,a,b}]

attempts to approximate ∫abf(x)dx.

Example 3.34

Evaluate ∫0π3e−x2cos⁡x3dx.

Solution

In this case, Mathematica is unable to evaluate the integral with Integrate.

We use the button to complete the Integrate command.

i1 = Integrate [ Exp [ − x ∧ 2 ] Cos [ x ∧ 3 ] , { x , 0 , Pi ∧ ( 1 / 3 ) } ]

∫ 0 π 1 / 3 e − x 2 Cos [ x 3 ] d x

An approximation is obtained with N.

N [ i1 ]

0.701566

Instead of using Integrate followed by N, you can use NIntegrate to numerically evaluate the integral.

NIntegrate [ Exp [ − x ∧ 2 ] Cos [ x ∧ 3 ] , { x , 0 , Pi ∧ ( 1 / 3 ) } ]

0.701566

returns the same result as that obtained using Integrate followed by N. □

In some cases, you may wish to investigate particular numerical methods that can be used to approximate integrals. To implement numerical methods like Simpson's rule or the trapezoidal rule, redefine the function leftsum (middlesum or rightsum) discussed previously to perform the calculation for the desired method.

3.3.4 Area

Suppose that y=f(x) and y=g(x) are continuous on [a,b] and that f(x)⩾g(x) for a⩽x⩽b. The area of the region bounded by the graphs of y=f(x), y=g(x), x=a, and x=b is

A = ∫ a b [ f ( x ) − g ( x ) ] d x .

(3.14)

Sometimes determining the “greater” function and the “lower” function can be difficult. Equation (3.14) in its more general form tells us that the region bounded by the graphs of y=f(x), y=g(x), x=a, and x=b is

A = ∫ a b | f ( x ) − g ( x ) | d x .

(3.15)

Example 3.35

Find the area between the graphs of y=sin⁡x and y=cos⁡x on the interval [0,2π].

Solution

We graph y=sin⁡x and y=cos⁡x on the interval [0,2π] in Fig. 3.36 with Plot. The graph of y=cos⁡x is dashed. Observe that including the option Filling->{1->{2}} fills the region between the two plots.

Figure 3.36 y=sin⁡x and y=cos⁡x on the interval [0,2π] (University of Wisconsin colors).

Plot [ { Sin [ x ] , Cos [ x ] } , { x , 0 , 2 π } ,

PlotStyle → { { CMYKColor [ . 03 , 1 , . 66 , . 12 ] , Thickness [ . 01 ] } ,

{ Thickness [ . 01 ] , Dashing [ { 0.025 } ] , CMYKColor [ . 05 , . 26 , 1 , . 27 ] } } , Filling → { 1 → { 2 } } ,

AspectRatio → Automatic ]

To find the upper and lower limits of integration, we must solve the equation sin⁡x=cos⁡x for x. Observe that Mathematica returns all solutions to the equation that it can find. The results show us that there are infinitely many solutions to the equation.

Solve [ Sin [ x ] == Cos [ x ] , x ]

{ { x → ConditionalExpression [ − 3 π 4 + 2 π C [ 1 ] , C [ 1 ] ∈ Integers ] } ,

{ x → ConditionalExpression [ π 4 + 2 π C [ 1 ] , C [ 1 ] ∈ Integers ] } }

For us the solutions of interest are valid for 0⩽x⩽2π, which are x=π/4 and x=5π/4. We check that these are valid solutions of sin⁡x=cos⁡x with ==; in each case the returned result is True.

Sin [ π 4 ] == Cos [ π 4 ]

Sin [ 5 π 4 ] == Cos [ 5 π 4 ]

True

True

Hence, the area of the region between the graphs is given by

A = ∫ 0 π / 4 [ cos ⁡ x − sin ⁡ x ] d x + ∫ π / 4 5 π / 4 [ sin ⁡ x − cos ⁡ x ] d x + ∫ 5 π / 4 2 π [ cos ⁡ x − sin ⁡ x ] d x .

(3.16)

Notice that if we take advantage of symmetry we can simplify (3.16) to

A = 2 ∫ π / 4 5 π / 4 [ sin ⁡ x − cos ⁡ x ] d x .

(3.17)

We evaluate (3.17) with Integrate to see that the area of the region between the two graphs is 42.

∫ 0 π 4 ( Cos [ x ] − Sin [ x ] ) d x + ∫ π 4 5 π 4 ( Sin [ x ] − Cos [ x ] ) d x + ∫ 5 π 4 2 π ( Cos [ x ] − Sin [ x ] ) d x

42  □

In cases when we cannot calculate the points of intersection of two graphs exactly, we can frequently use FindRoot to approximate the points of intersection.

Example 3.36

Let

p ( x ) = 3 10 x 5 − 3 x 4 + 11 x 3 − 18 x 2 + 12 x + 1

and

q ( x ) = − 4 x 3 + 28 x 2 − 56 x + 32 .

Approximate the area of the region bounded by the graphs of y=p(x) and y=q(x).

Solution

After defining p and q, we graph them on the interval [−1,5] in Fig. 3.37 to obtain an initial guess of the intersection points of the two graphs.

Figure 3.37 p and q on the interval [−1,5] (University of Maryland colors).

When you use Tooltip, you can slide your cursor over a plot and the function being graphed is displayed.

Clear [ p , q ]

p [ x _ ] = 3 x 5 10 − 3 x 4 + 11 x 3 − 18 x 2 + 12 x + 1 ;

q [ x _ ] = − 4 x 3 + 28 x 2 − 56 x + 32 ;

Plot [ Tooltip [ { p [ x ] , q [ x ] } ] , { x , − 1 , 5 } ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 0 , . 91 , . 76 , . 06 ] } ,

{ Thickness [ . 01 ] , CMYKColor [ 0 , . 15 , . 94 , 0 ] } } ,

AxesLabel → { x , y } , AspectRatio → 1 ]

The x-coordinates of the three intersection points are the solutions of the equation p(x)=q(x). Although Mathematica can solve this equation exactly, approximate solutions are more useful for the problem and obtained with NSolve.

intpts = NRoots [ p [ x ] == q [ x ] , x ]

x = = 0.772058 ‖ x = = 1.5355 − 3.57094 i ‖ x = = 1.5355 + 3.57094 i ‖ x = = 2.29182 ‖ x = = 3.86513

The numbers are extracted from the list with Part ([[...]]). For example, 0.772058 is the second part of the first part of intpts. Counting from left to right, 2.29182 is the second part of the fourth part of intpts.

x1 = intpts [ [ 1 , 2 ] ]

x2 = intpts [ [ 4 , 2 ] ]

x3 = intpts [ [ 5 , 2 ] ]

0.772058

2.29182

3.86513

Using the roots to the equation p(x)=q(x) and the graph we see that p(x)⩾q(x) for 0.772⩽x⩽2.292 and q(x)⩾p(x) for 2.292⩽x⩽3.865. Hence, an approximation of the area bounded by p(x) and q(x) is given by the sum

∫ 0.772 2.292 [ p ( x ) − q ( x ) ] d x + ∫ 2.292 3.865 [ q ( x ) − p ( x ) ] d x .

These two integrals are computed with Integrate and NIntegrate. As expected, the two values are the same.

∫ x1 x2 ( p [ x ] − q [ x ] ) d x + ∫ x2 x3 ( q [ x ] − p [ x ] ) d x

12.1951

NIntegrate [ p [ x ] − q [ x ] , { x , x1 , x2 } ] + NIntegrate [ q [ x ] − p [ x ] , { x , x2 , x3 } ]

12.1951

We conclude that the area is approximately 12.195 units². □

Parametric Equations

If the curve, C, defined parametrically by x=x(t), y=y(t), a⩽t⩽b is a nonnegative continuous function of x and x(a)<x(b) the area under the graph of C and above the x-axis is

∫ x ( a ) x ( b ) y d x = ∫ a b y ( t ) x ′ ( t ) d t .

Graphically, y is a function of x, y=y(x), if the graph of y=y(x) passes the vertical line test.

Example 3.37

The Astroid

Find the area enclosed by the astroid x=sin3⁡t, y=cos3⁡t, 0⩽t⩽2π.

Solution

We begin by defining x and y and then graphing the astroid with ParametricPlot in Fig. 3.38.

Figure 3.38 The astroid x=sin3⁡t, y=cos3⁡t, 0 ⩽ t ⩽ 2π (Harvard University colors).

x [ t _ ] = Sin [ t ] ∧ 3 ;

y [ t _ ] = Cos [ t ] ∧ 3 ;

ParametricPlot [ { x [ t ] , y [ t ] } , { t , 0 , 2 Pi } , AspectRatio - > Automatic ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ . 07 , . 94 , . 65 , . 25 ] } ]

Observe that x(0)=0 and x(π/2)=1 and the graph of the astroid in the first quadrant is given by x=sin3⁡t, y=cos3⁡t, 0⩽t⩽π/2. Hence, the area of the astroid in the first quadrant is given by

∫ 0 π / 2 y ( t ) x ′ ( t ) d t = 3 ∫ 0 π / 2 sin 2 ⁡ t cos 4 ⁡ t d t

and the total area is given by

A = 4 ∫ 0 π / 2 y ( t ) x ′ ( t ) d t = 12 ∫ 0 π / 2 sin 2 ⁡ t cos 4 ⁡ t d t = 3 8 π ≈ 1.178 ,

which is computed with Integrate and then approximated with N.

area = 4 Integrate [ y [ t ] x ′ [ t ] , { t , 0 , Pi / 2 } ]

3 π 8

N [ area ]

1.1781  □

Polar Coordinates

For problems involving “circular symmetry” it is often easier to work in polar coordinates. The relationship between (x,y) in rectangular coordinates and (r,θ) in polar coordinates is given by

x = r cos ⁡ θ y = r sin ⁡ θ

and

r 2 = x 2 + y 2 tan ⁡ θ = y x .

If r=f(θ) is continuous and nonnegative for α⩽θ⩽β, then the area A of the region enclosed by the graphs of r=f(θ), θ=α, and θ=β is

A = 1 2 ∫ α β [ f ( θ ) ] 2 d θ = 1 2 ∫ α β r 2 d θ .

Example 3.38

Lemniscate of Bernoulli

The lemniscate of Bernoulli is given by

( x 2 + y 2 ) 2 = a 2 ( x 2 − y 2 ) ,

where a is a constant. (a) Graph the lemniscate of Bernoulli if a=2. (b) Find the area of the region bounded by the lemniscate of Bernoulli.

Solution

This problem is much easier solved in polar coordinates so we first convert the equation from rectangular to polar coordinates with ReplaceAll (/.) and then solve for r with Solve.

lofb = ( x ∧ 2 + y ∧ 2 ) ∧ 2 == a ∧ 2 ( x ∧ 2 − y ∧ 2 ) ;

topolar=lofb/.{x->rCos[t],y->rSin[t]}

( r 2 Cos [ t ] 2 + r 2 Sin [ t ] 2 ) 2 = = a 2 ( r 2 Cos [ t ] 2 − r 2 Sin [ t ] 2 )

Solve [ topolar , r ] // Simplify

{ { r → 0 } , { r → 0 } , { r → − a 2 Cos [ 2 t ] } , { r → a 2 Cos [ 2 t ] } }

These results indicate that an equation of the lemniscate in polar coordinates is r2=a2cos⁡2θ. The graph of the lemniscate is then generated in Fig. 3.39 (top figure) using PolarPlot. The portion of the lemniscate in quadrant one is obtained by graphing r=2cos⁡2θ, 0⩽θ⩽π/4.

Figure 3.39 On the left, the lemniscate and on the portion of the lemniscate in quadrant 1 (Harvard University colors).

p1 = PolarPlot [ { − 2 Sqrt [ Cos [ 2 t ] ] , 2 Sqrt [ Cos [ 2 t ] ] } , { t , 0 , 2 Pi } ,

PlotStyle → { Thickness [ . 015 ] , CMYKColor [ . 16 , . 23 , . 23 , . 44 ] } ,

AxesLabel → { x , y } ] ;

p2 = PolarPlot [ 2 Sqrt [ Cos [ 2 t ] ] , { t , 0 , Pi / 4 } ,

PlotStyle → { Thickness [ . 015 ] , CMYKColor [ . 33 , . 15 , . 14 , . 17 ] } ,

AxesLabel → { x , y } ] ;

Show [ GraphicsRow [ { p1 , p2 } ] ]

Then, taking advantage of symmetry, the area of the lemniscate is given by

A = 2 ⋅ 1 2 ∫ − π / 4 π / 4 r 2 d θ = 2 ∫ 0 π / 4 r 2 d θ = 2 ∫ 0 π / 4 a 2 cos ⁡ 2 θ d θ = a 2 ,

which we calculate with Integrate.

Integrate [ 2 a ∧ 2 Cos [ 2 t ] , { t , 0 , Pi / 4 } ]

a2  □

3.3.5 Arc Length

Let y=f(x) be a function for which f′(x) is continuous on an interval [a,b]. Then the arc length of the graph of y=f(x) from x=a to x=b is given by

L = ∫ a b ( d y d x ) 2 + 1 d x .

(3.18)

The resulting definite integrals used for determining arc length are usually difficult to compute because they involve a radical. In these situations, Mathematica is helpful with approximating solutions to these types of problems.

Example 3.39

Find the length of the graph of y=x48+14x2 from (a) x=1 to x=2 and from (b) x=−2 to x=−1.

Solution

With no restrictions on the value of x, x2=|x|. Generally, Mathematica does not automatically algebraically simplify (dy/dx)2+1 because Mathematica does not know if x is positive or negative.

y [ x _ ] = x ∧ 4 / 8 + 1 / ( 4 x ∧ 2 ) ;

i1 = Factor [ y ′ [ x ] ∧ 2 + 1 ]

( 1 + x 2 ) 2 ( 1 − x 2 + x 4 ) 2 4 x 6

i2 = PowerExpand [ Sqrt [ i1 ] ]

(1+x2)(1−x2+x4)2x3

PowerExpand[expr] simplifies radicals in the expression expr assuming that all variables are positive.

In fact, for (b), x is negative so 12(x6+1)2x6=−12x6+1x3. Mathematica simplifies 12(x6+1)2x6=12x6+1x3 and correctly evaluates the arc length integral (3.18) for (a).

Integrate [ Sqrt [ y ′ [ x ] ∧ 2 + 1 ] , { x , 1 , 2 } ]

33 16

For (b), we compute the arc length integral (3.18).

Integrate [ Sqrt [ y ′ [ x ] ∧ 2 + 1 ] , { x , − 2 , − 1 } ]

33 16

As we expect, both values are the same. □

Parametric Equations

If the smooth curve, C, defined parametrically by x=x(t), y=y(t), t∈[a,b] is traversed exactly once as t increases from t=a to t=b, the arc length of C is given by

L = ∫ a b ( d x d t ) 2 + ( d y d t ) 2 d t .

(3.19)

C is smooth if both x′(t) and y′(t) are continuous on (a,b) and not simultaneously zero for t∈(a,b).

Example 3.40

Find the length of the graph of x=2t2, y=2t−12t3, −2⩽t⩽2.

Solution

For illustrative purposes, we graph x=2t2, y=2t−12t3 for −3⩽t⩽3 and −2⩽t⩽2 (thickened) in Fig. 3.40.

Figure 3.40 x=2t2, y=2t−12t3 (Pennsylvania State University colors).

x [ t _ ] = t ∧ 2 Sqrt [ 2 ] ; y [ t _ ] = 2 t − 1 / 2 t ∧ 3 ;

p1 = ParametricPlot [ { x [ t ] , y [ t ] } , { t , − 3 , 3 } ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ . 75 , . 2 , 0 , 0 ] } ] ;

p2 = ParametricPlot [ { x [ t ] , y [ t ] } , { t , − 2 , 2 } ,

PlotStyle → { CMYKColor [ . 75 , . 66 , 0 , 0 ] , Thickness [ . 02 ] } ] ;

Show [ p1 , p2 , PlotRange - > All , AxesLabel → { x , y } ]

Mathematica is able to compute the exact value of the arc length (3.19) although the result is quite complicated. For length considerations, the result of entering the i1 command are not displayed here.

Factor [ x ′ [ t ] ∧ 2 + y ′ [ t ] ∧ 2 ]

14(4−4t+3t2)(4+4t+3t2)

i1 = Integrate [ 2 Sqrt [ x ′ [ t ] ∧ 2 + y ′ [ t ] ∧ 2 ] , { t , 0 , 2 } ]

∫ 0 2 2 8 t 2 + ( 2 − 3 t 2 2 ) 2 d t

A more meaningful approximation is obtained with N or using NIntegrate.

N [ i1 ]

13.7099

NIntegrate [ 2 Sqrt [ x ′ [ t ] ∧ 2 + y ′ [ t ] ∧ 2 ] , { t , 0 , 2 } ]

13.7099

We conclude that the arc length is approximately 13.71.

Observe that Mathematica 11 cannot evaluate the definite integral directly. However, if we first compute an anti-derivative with integrate in i2,

i2 = Integrate [ 2 Sqrt [ x ′ [ t ] ∧ 2 + y ′ [ t ] ∧ 2 ] , t ]

( 27 i i − 2 2 t ( 16 + 8 t 2 + 9 t 4 ) − 16 ( i + 2 2 ) 4 i − 8 2 + 9 i t 2 i − 2 2 4 i + 8 2 + 9 i t 2 i + 2 2

EllipticE [ i ArcSinh [ 3 2 i i − 2 2 t ] , i − 2 2 i + 2 2 ] + 32 ( − 4 i + 2 ) 4 i − 8 2 + 9 i t 2 i − 2 2 4 i + 8 2 + 9 i t 2 i + 2 2

EllipticF [ i ArcSinh [ 3 2 i i − 2 2 t ] , i − 2 2 i + 2 2 ] ) / ( 81 i i − 2 2 16 + 8 t 2 + 9 t 4 )

and then apply the Fundamental Theorem of Calculus by subtracting the value of i2 if t=0 from the value of i2 if t=2,

ul = i2 /. t → 2

ll = i2 /. t → 0

val = ul − ll

1 648 3 i i − 2 2 ( 10368 i i − 2 2 − 16 40 i − 8 2 i − 2 2 ( i + 2 2 ) 40 i + 8 2 i + 2 2

EllipticE [ i ArcSinh [ 3 i i − 2 2 ] , i − 2 2 i + 2 2 ] + 32 40 i − 8 2 i − 2 2 ( − 4 i + 2 ) 40 i + 8 2 i + 2 2

EllipticF [ i ArcSinh [ 3 i i − 2 2 ] , i − 2 2 i + 2 2 ] )

0

1 648 3 i i − 2 2 ( 10368 i i − 2 2 − 16 40 i − 8 2 i − 2 2 ( i + 2 2 ) 40 i + 8 2 i + 2 2

EllipticE [ i ArcSinh [ 3 i i − 2 2 ] , i − 2 2 i + 2 2 ] + 32 40 i − 8 2 i − 2 2 ( − 4 i + 2 ) 40 i + 8 2 i + 2 2

EllipticF[iArcSinh[3ii−22],i−22i+22])

N [ val ]

13.7099 + 1.7763568394002505 ` * − ∧ 15 i

N [ val ] // Chop

13.7099

we obtain the same result. □

Polar Coordinates

If the smooth polar curve C given by r=f(θ), α⩽θ⩽β is traversed exactly once as θ increases from α to β, the arc length of C is given by

L = ∫ α β ( d r d θ ) 2 + r 2 d θ

(3.20)

Example 3.41

Find the length of the graph of r=θ, 0⩽θ⩽10π.

Solution

We begin by defining r and then graphing r with PolarPlot in Fig. 3.41.

Figure 3.41 r = θ for 0 ⩽ θ ⩽ 10π (Pennsylvania State University colors).

r [ t _ ] = t ;

PolarPlot [ r [ t ] , { t , 0 , 10 Pi } , AspectRatio - > Automatic ,

PlotStyle → { CMYKColor [ 0 , . 75 , . 6 , 0 ] , Thickness [ . 01 ] } ,

AxesLabel→{x,y}]

Using (3.20), the length of the graph of r is given by ∫010π1+θ2dθ. The exact value is computed with Integrate

ev = Integrate [ Sqrt [ r ′ [ t ] ∧ 2 + r [ t ] ∧ 2 ] , { t , 0 , 10 Pi } ]

5 π 1 + 100 π 2 + 1 2 ArcSinh [ 10 π ]

and then approximated with N.

N [ ev ]

495.801

We conclude that the length of the graph is approximately 495.8 units. □

3.3.6 Solids of Revolution

Volume

Let y=f(x) be a nonnegative continuous function on [a,b]. The volume of the solid of revolution obtained by revolving the region bounded by the graphs of y=f(x), x=a, x=b, and the x-axis about the x-axis is given by

V = π ∫ a b [ f ( x ) ] 2 d x .

(3.21)

If 0⩽a<b, the volume of the solid of revolution obtained by revolving the region bounded by the graphs of y=f(x), x=a, x=b, and the x-axis about the y-axis is given by

V = 2 π ∫ a b x f ( x ) d x .

(3.22)

Example 3.42

Let g(x)=xsin2⁡x. Find the volume of the solid obtained by revolving the region bounded by the graphs of y=g(x), x=0, x=π, and the x-axis about (a) the x-axis; and (b) the y-axis.

Solution

After defining g, we graph g on the interval [0,π] in Fig. 3.42 (a).

Figure 3.42 (a) g(x) for 0 ⩽ x ⩽ π. (b) g(x) revolved about the x-axis. (c) g(x) revolved about the y-axis (University of California–Santa Barbara colors).

With Mathematica 11, for three dimensional graphics, you can adjust the viewpoint by clicking on the three-dimensional graphics object and dragging to the desired viewing angle.

g [ x _ ] = x Sin [ x ] ∧ 2 ;

p1 = Plot [ g [ x ] , { x , 0 , Pi } , AspectRatio - > Automatic ,

PlotLabel → “(a)” ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 1 , . 82 , . 17 , . 04 ] } } ] ;

The volume of the solid obtained by revolving the region about the x-axis is given by equation (3.21) while the volume of the solid obtained by revolving the region about the y-axis is given by equation (3.22). These integrals are computed with Integrate and named xvol and yvol, respectively. We use N to approximate each volume.

xvol = Integrate [ Pi g [ x ] ∧ 2 , { x , 0 , Pi } ]

N [ xvol ]

1 64 π 2 ( − 15 + 8 π 2 )

9.86295

yvol = Integrate [ 2 Pi x g [ x ] , { x , 0 , Pi } ]

N [ yvol ]

1 6 π 2 ( − 3 + 2 π 2 )

27.5349

We can use ParametricPlot3D to visualize the resulting solids by parametrically graphing the equations given by

{ x = r cos ⁡ t y = r sin ⁡ t z = g ( r )

for r between 0 and π and t between −π and π to visualize the graph of the solid obtained by revolving the region about the y-axis and by parametrically graphing the equations given by

{ x = r y = g ( r ) cos ⁡ t z = g ( r ) sin ⁡ t

for r between 0 and π and t between −π and π to visualize the graph of the solid obtained by revolving the region about the x-axis. (See Figs. 3.42 (b) and 3.42 (c).) In this case, we identify the z-axis as the y-axis. Notice that we are simply using polar coordinates for the x and y-coordinates, and the height above the x, y-plane is given by z=g(r) because r is replacing x in the new coordinate system.

p2 = ParametricPlot3D [ { r , g [ r ] Cos [ t ] , g [ r ] Sin [ t ] } , { r , 0 , Pi } , { t , 0 , 2 Pi } ,

PlotPoints - > { 30 , 30 } , PlotLabel → “(b)” ,

PlotStyle → { CMYKColor [ 0 , . 18 , . 89 , 0 ] , Specularity [ White , 10 ] } ,

Mesh→False];

p3 = ParametricPlot3D [ { r Cos [ t ] , r Sin [ t ] , g [ r ] } , { r , 0 , Pi } , { t , 0 , 2 Pi } ,

PlotPoints - > { 30 , 30 } , PlotLabel → “(c)” ,

PlotStyle → { CMYKColor [ 1 , . 86 , . 36 , . 28 ] ,

Opacity [ . 9 ] , Specularity [ White , 1 ] } ,

Mesh → False ] ;

p1, p2, and p3 are shown together side-by-side in Fig. 3.42 using Show together with GraphicsRow.

Show[GraphicsRow[{p1,p2,p3}]]  □

We now demonstrate a volume problem that requires the method of disks.

Example 3.43

Let f(x)=e−(x−3)2cos⁡[4(x−3)]. Approximate the volume of the solid obtained by revolving the region bounded by the graphs of y=f(x), x=1, x=5, and the x-axis about the x-axis.

Solution

Proceeding as in the previous example, we first define and graph f on the interval [1,5] in Fig. 3.43 (a).

Figure 3.43 (a) f(x) for 1 ⩽ x ⩽ 5. (b) f(x) revolved abut the x-axis (University of Southern California colors).

f [ x _ ] = Exp [ − ( x − 3 ) ∧ 2 Cos [ 4 ( x − 3 ) ] ] ;

p1 = Plot [ f [ x ] , { x , 1 , 5 } , AspectRatio - > Automatic ,

PlotLabel → “(a)” , AxesLabel → { x , y } ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ . 07 , 1 , . 65 , . 32 ] } ] ;

In this case, an approximation is desired so we use NIntegrate to approximate the integral V=∫15π[f(x)]2dx.

NIntegrate[Pif[x]∧2,{x,1,5}]

16.0762

In the same manner as before, ParametricPlot3D can be used to visualize the resulting solid by graphing the set of equations given parametrically by

{ x = r y = f ( r ) cos ⁡ t z = f ( r ) sin ⁡ t

for r between 1 and 5 and t between 0 and 2π. In this case, polar coordinates are used in the y, z-plane with the distance from the x-axis given by f(x). Because r replaces x in the new coordinate system, f(x) becomes f(r) in these equations. See Fig. 3.43 (b).

p2 = ParametricPlot3D [ { r , f [ r ] Cos [ t ] , f [ r ] Sin [ t ] } , { r , 1 , 5 } , { t , 0 , 2 Pi } ,

PlotPoints - > { 45 , 35 } , PlotLabel → “(b)” ,

PlotStyle → { CMYKColor [ 0 , . 27 , 1 , 0 ] } ] ;

Show[GraphicsRow[{p1,p2}]]  □

When revolving a curve about the y-axis, you can use RevolutionPlot3D rather than the parametrization given previously.

Example 3.44

Let f(x)=exp⁡(−2(x−2)2)+exp⁡(−(x−4)2) for 0⩽x⩽6. (a) Find the minimum and maximum values of f(x) on [0,6]. Let R be the region bounded by y=f(x), x=0, x=6, and the y-axis. (b) Find the volume of the solid obtained by revolving R about the y-axis. (c) Find the volume of the solid obtained by revolving R about the x-axis.

Solution

(a) Although Maximize and Minimize cannot find the exact maximum and minimum values, using N or NMaximize and NMinimize give accurate approximations.

NMaximize and NMinimize work in the same way as Maximize and Minimize but return approximations rather than exact results.

f [ x _ ] = Exp [ − 2 ( x − 2 ) ∧ 2 ] + Exp [ − ( x − 4 ) ∧ 2 ] ;

Maximize [ f [ x ] , x ]

{ e − ( − 4 + Root [ { − 4 + 2 # 1 + e # 1 2 ( − 4 e 8 + # 1 e 8 ) & , 2.0196244769513376905 } ] ) 2

× e − 2 ( − 2 + Root [ { − 4 + 2 # 1 + e # 1 2 ( − 4 e 8 + # 1 e 8 ) & , 2.0196244769513376905 } ] ) 2

( e ( − 4 + Root [ { − 4 + 2 # 1 + e # 1 2 ( − 4 e 8 + # 1 e 8 ) & , 2.0196244769513376905 } ] ) 2

+ e 2 ( − 2 + Root [ { − 4 + 2 # 1 + e # 1 2 ( − 4 e 8 + # 1 e 8 ) & , 2.0196244769513376905 } ] ) 2 ) , { x → Root [ { − 4 + 2 # 1

+ e # 1 2 ( − 4 e 8 + # 1 e 8 ) & , 2.0196244769513376905 } ] } }

Maximize [ f [ x ] , x ] // N

{1.01903,{x→2.01962}}

NMaximize [ { f [ x ] , 0 ⩽ x ⩽ 6 } , x ]

{ 1.00034 , { x → 3.99864 } }

Minimize [ { f [ x ] , 0 ⩽ x ⩽ 6 } , x ]

{ e 8 + e 16 e 24 , { x → 0 } }

Minimize [ { f [ x ] , 0 ⩽ x ⩽ 6 } , x ] // N

{ 0.000335575 , { x → 0 . } }

NMinimize [ { f [ x ] , 0 ⩽ x ⩽ 6 } , x ] // N

{ 0.000335575 , { x → 0 . } }

We double check these results by graphing f(x) and f′(x) in Fig. 3.44 and then using FindRoot to approximate the critical numbers.

Figure 3.44 We use the graph of f^′(x) to help us estimate the initial values to approximate the critical numbers with FindRoot (University of Minnesota colors).

pf1 = Plot [ f [ x ] , { x , 0 , 6 } , PlotLabel → “(a)” ,

AxesLabel → { x , y } , PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 0 , . 16 , 1 , 0 ] } } ]

pf2 = Plot [ Tooltip [ { f [ x ] , f ′ [ x ] } ] , { x , 0 , 6 } ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ 0 , . 16 , 1 , 0 ] } ,

{ Thickness [ . 01 ] , CMYKColor [ 0 , 1 , . 63 , . 29 ] } } ,

PlotLabel → “(b)” , AxesLabel → { x , y } ]

Show [ GraphicsRow [ { pf1 , pf2 } ] ]

Map [ FindRoot [ f ′ [ x ] = = 0 , { x , # 1 } ] & , { 2 , 3 , 4 } ]

{{x→2.01962},{x→2.92167},{x→3.99864}}

(b) Mathematica finds the exact volume of the solids although the results are expressed in terms of the Error function, Erf.

Integrate [ Pi x f [ x ] , { x , 0 , 6 } ]

π ( − 1 + 2 e 16 + e 24 − 2 e 28 + 2 e 32 π ( 4 Erf [ 2 ] + 4 Erf [ 4 ] + 2 ( Erf [ 2 2 ] + Erf [ 4 2 ] ) ) ) 4 e 32

NIntegrate [ Pi x f [ x ] , { x , 0 , 6 } ]

30.0673

Integrate [ Pi f [ x ] ∧ 2 , { x , 0 , 6 } ]

1 12 e 8 / 3 π 3 / 2 ( 3 e 8 / 3 ( Erf [ 4 ] + Erf [ 8 ] + 2 ( Erf [ 2 2 ] + Erf [ 4 2 ] ) ) + 4 3 ( Erf [ 8 3 ] + Erf [ 10 3 ] ) )

NIntegrate [ Pi f [ x ] ∧ 2 , { x , 0 , 6 } ]

7.1682

To visualize the solid revolved about the y-axis, we use RevolutionPlot3D in p1. We generate the curve in p2, a set of axes, and a representative “slice” of the curve. See Fig. 3.45 (a). After, we show the solid together with a representative shell. See Fig. 3.45 (b).

Figure 3.45 (a) The solid. (b) The solid with a “typical” shell. (c) Several shells.

p1 = RevolutionPlot3D [ f [ x ] , { x , 1 , 5 } ,

BoxRatios → { 2 , 2 , 1 } , PlotRange → { { − 5 , 5 } , { − 5 , 5 } , { 0 , 5 / 4 } } ,

Mesh → None , PlotStyle → Opacity [ . 4 ] ,

ColorFunction → “LightTemperatureMap” ] ;

p2 = ParametricPlot3D [ { x , 0 , Exp [ − 2 ( x − 2 ) ∧ 2 ] + Exp [ − ( x − 4 ) ∧ 2 ] } ,

{ x , 1 , 5 } , { t , 0 , 2 Pi } ,

PlotStyle → Thickness [ . 05 ] , BoxRatios → { 2 , 2 , 1 } ,

Axes → Automatic , Boxed → False ] ;

p3 = ParametricPlot3D [ { x , 0 , 0 } , { x , − 5 , 5 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 075 ] } , BoxRatios → { 2 , 2 , 1 } ,

Axes → Automatic , Boxed → False ] ;

p4 = ParametricPlot3D [ { 0 , 0 , x } , { x , 0 , 5 / 4 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 1 ] } , BoxRatios → { 2 , 2 , 1 } ,

Axes → Automatic , Boxed → False ] ;

p5 = Graphics3D [ { Gray , Thickness [ . 01 ] , Line [ { { 3.6 , 0 , 0 } ,

{ 3.6 , 0 , Exp [ − 2 ( 3.6 − 2 ) ∧ 2 ] + Exp [ − ( 3.6 − 4 ) ∧ 2 ] } } ] } ] ;

p6 = ParametricPlot3D [ { 3.6 Cos [ t ] , 3.6 Sin [ t ] , z } , { t , 0 , 2 Pi } ,

{ z , 0 , Exp [ − 2 ( 3.6 − 2 ) ∧ 2 ] + Exp [ − ( 3.6 − 4 ) ∧ 2 ] } , Mesh → None ,

PlotStyle → Opacity [ . 8 ] , ColorFunction → “TemperatureMap” ] ;

g1 = Show [ p1 , p2 , p3 , p4 , p5 , Boxed → False , Axes → None ]

g2 = Show [ p1 , p2 , p3 , p4 , p5 , p6 , Boxed → False , Axes → None ]

Show [ GraphicsRow [ { g1 , g2 } ] ]

Finally, we show the solid together with several shells in Fig. 3.45 (c).

sp = Table [ j // N , { j , 1.1 , 4.9 , 3.8 / 14 } ] ;

p7 = Table [ ParametricPlot3D [ { sp [ [ i ] ] Cos [ t ] , sp [ [ i ] ] Sin [ t ] , z } ,

{ t , 0 , 2 Pi } ,

{ z , 0 , Exp [ − 2 ( sp [ [ i ] ] − 2 ) ∧ 2 ] + Exp [ − ( sp [ [ i ] ] − 4 ) ∧ 2 ] } ,

Mesh → None ,

PlotStyle → Opacity [ . 8 ] ,

ColorFunction → “TemperatureMap” ] , { i , 1 , Length [ sp ] } ] ;

g3 = Show [ p1 , p2 , p3 , p4 , p5 , p7 , Boxed → False , Axes → None ]

Show [ GraphicsRow [ { g1 , g2 , g3 } ] ]

You can use Animate to animate the process as illustrated in Fig. 3.46.

Figure 3.46 Using Animate to illustrate the finding the volume of a solid using the method of shells.

Animate [

p8 = ParametricPlot3D [ { i Cos [ t ] , i Sin [ t ] , z } , { t , 0 , 2 Pi } ,

{ z , 0 , Exp [ − 2 ( i − 2 ) ∧ 2 ] + Exp [ − ( i − 4 ) ∧ 2 ] } , Mesh → None ,

PlotStyle → Opacity [ . 8 ] , ColorFunction → “TemperatureMap” ] ;

Show [ p1 , p2 , p3 , p4 , p5 , p8 , Boxed → False , Axes → None ] ,

{i,1.1,4.9}]

For revolving f(x) about the x-axis, we proceed in much the same way. First, we plot f(x) with a set of axes in three-space.

f [ x _ ] = Exp [ − 2 ( x − 2 ) ∧ 2 ] + Exp [ − ( x − 4 ) ∧ 2 ] ;

p1 = ParametricPlot3D [ { x , 0 , f [ x ] } , { x , 0 , 6 } , PlotStyle → { Thick , Black } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

p1b = ParametricPlot3D [ { x , 0 , − f [ x ] } , { x , 0 , 6 } , PlotStyle → { Thick , Black } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

p2 = ParametricPlot3D [ { x , 0 , 0 } , { x , 0 , 6 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 075 ] } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

p3 = ParametricPlot3D [ { 0 , 0 , x } , { x , − 3 / 2 , 3 / 2 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 1 ] } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

Show [ p1 , p1b , p2 , p3 ]

(c) Next, we generate a basic plot of the solid in p4 and then a set of disks inside the solid in t3d.

f [ x _ ] = Exp [ − 2 ( x − 2 ) ∧ 2 ] + Exp [ − ( x − 4 ) ∧ 2 ] ;

p1 = ParametricPlot3D [ { x , 0 , f [ x ] } , { x , 0 , 6 } , PlotStyle → { Thick , Black } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

p1b = ParametricPlot3D [ { x , 0 , − f [ x ] } , { x , 0 , 6 } , PlotStyle → { Thick , Black } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } , BoxRatios → { 1 , 1 , 1 } ] ;

p2 = ParametricPlot3D [ { x , 0 , 0 } , { x , 0 , 6 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 075 ] } , PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } ] ;

p3 = ParametricPlot3D [ { 0 , 0 , x } , { x , − 3 / 2 , 3 / 2 } , { t , 0 , 2 Pi } ,

PlotStyle → { Gray , Thickness [ . 1 ] } , PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } ] ;

Show [ p1 , p1b , p2 , p3 ]

p4 = ParametricPlot3D [ { r , f [ r ] Cos [ t ] , f [ r ] Sin [ t ] } ,

{ r , 0 , 6 } , { t , 0 , 2 Pi } , PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } ]

t3d = Table [ ParametricPlot3D [ { x , r f [ x ] Cos [ t ] , r f [ x ] Sin [ t ] } , { r , 0 , 1 } , { t , 0 , 2 Pi } ,

PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } , ColorFunction → “TemperatureMap” , Mesh → 5 ] , { x , 0 , 6 , 6 / 14 } ] ;

Two variations of the solid are plotted in p5 and p6. In each case, we use MeshFunctions to have the contour lines (mesh) correspond to f(x) values rather than the rectangular default mesh. In p6 the solid is made transparent with the Opacity option.

p5 = ParametricPlot3D [ { r , f [ r ] Cos [ t ] , f [ r ] Sin [ t ] } ,

{ r , 0 , 6 } , { t , 0 , 2 Pi } , PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } , { − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } , MeshFunctions - > { # 1 & } , Mesh → 60 ]

p6 = ParametricPlot3D [ { r , f [ r ] Cos [ t ] , f [ r ] Sin [ t ] } ,

{ r , 0 , 6 } , { t , 0 , 2 Pi } , PlotRange → { { 0 , 6 } , { − 3 / 2 , 3 / 2 } ,

{ − 3 / 2 , 3 / 2 } } ,

BoxRatios → { 1 , 1 , 1 } , MeshFunctions - > { # 1 & } , Mesh → 25 ,

PlotStyle → Opacity [ . 2 ] ,

MeshStyle → { Gray , Thick } ] ;

Show [ p1 , p1b , p2 , p3 , p6 ]

Several combinations of the images are shown in Figs. 3.47 and 3.48.

Figure 3.47 (a) Seeing f(x) on the solid. (b) Disks in the solid.

Figure 3.48 (a) f(x) in space. (b) The basic solid. (c) Contours based on f(x) values. (d) Seeing f(x) on the solid.

Show [ GraphicsRow [ { Show [ p1 , p1b , p2 , p3 , p6 ] , Show [ p1 , p1b , p2 , p3 , p6 , t3d ] } ] ]

Show [ GraphicsGrid [ { { Show [ p1 , p1b , p2 , p3 ] , p4 } ,

{p5,Show[p1,p1b,p2,p3,p6]}}]]  □

To help identify regions, RegionPlot[constraints,{x,a,b},{y,a,b}] attempts to shade the region in the rectangle [a,b]×[c,d] that satisfies the constraints in constraints.

Example 3.45

Let g(x)=x, h(x)=x2, and R the region bounded by the graphs of g(x) and h(x). Find the volume of the solid obtained by revolving R about (a) the x-axis and (b) the y-axis.

Solution

We illustrate the use of RegionPlot to help us see R. See Fig. 3.49

Figure 3.49 (a) Graphs of f(x) and g(x). (b) The region in [0,2]×[0,2] for which x2⩽y⩽x. (c) The two plots displayed together.

g [ x _ ] = Sqrt [ x ] ;

h [ x _ ] = x ∧ 2 ;

p1a = Plot [ Tooltip [ { g [ x ] , h [ x ] } ] , { x , 0 , 2 } , PlotRange → { { 0 , 2 } , { 0 , 2 } } ,

AspectRatio → Automatic , PlotLabel → “(a)” ,

PlotStyle → { { CMYKColor [ . 85 , 1 , 0 , . 15 ] } , { CMYKColor [ . 5 , . 5 , . 5 , 1 ] } } ] ;

p1b = RegionPlot [ x ∧ 2 ⩽ y ⩽ Sqrt [ x ] , { x , 0 , 2 } , { y , 0 , 2 } , PlotStyle → Purple ,

PlotLabel → “(b)” ] ;

Show [ GraphicsRow [ { p1a , p1b , Show [ { p1a , p1b } , PlotLabel → “(c)” ] } ] ]

We plot the solids with ParametricPlot3D and contour lines along the function values using the MeshFunctions option in Fig. 3.50.

Figure 3.50 (a) The solid formed by revolving R about the x-axis. (b) The solid formed by revolving R about the y-axis (Northwestern University colors).

p4 = ParametricPlot3D [ { { r , g [ r ] Cos [ t ] , g [ r ] Sin [ t ] } ,

{ r , h [ r ] Cos [ t ] , h [ r ] Sin [ t ] } } , { r , 0 , 1 } ,

{ t , 0 , 2 Pi } , PlotRange → { { 0 , 3 / 2 } , { − 5 / 4 , 5 / 4 } , { − 5 / 4 , 5 / 4 } } ,

BoxRatios → { 1 , 1 , 1 } , MeshFunctions → { # 1 & } , PlotStyle → { Purple , Opacity [ . 1 ] } ,

PlotLabel → “(a)” ] ;

p5 = ParametricPlot3D [ { { r Cos [ t ] , r Sin [ t ] , g [ r ] } ,

{ r Cos [ t ] , r Sin [ t ] , h [ r ] } } , { r , 0 , 1 } ,

{ t , 0 , 2 Pi } , PlotRange → { { − 5 / 4 , 5 / 4 } , { − 5 / 4 , 5 / 4 } , { − 1 / 4 , 9 / 4 } } ,

BoxRatios → { 1 , 1 , 1 } , MeshFunctions → { # 1 & } , PlotStyle → { Purple , Opacity [ . 1 ] } ,

PlotLabel → “(b)” ] ;

Show [ GraphicsRow [ { p4 , p5 } ] ]

The volume of each solid is then found with Integrate and approximated with N.

Integrate [ Pi ( g [ x ] ∧ 2 − h [ x ] ∧ 2 ) , { x , 0 , 1 } ]

3 π 10

N [ % ]

0.942478

Integrate [ Pi x ( g [ x ] − h [ x ] ) , { x , 0 , 1 } ]

3 π 20

N [ % ]

0.471239 □

Surface Area

Let y=f(x) be a nonnegative function for which f′(x) is continuous on an interval [a,b]. Then the surface area of the solid of revolution obtained by revolving the region bounded by the graphs of y=f(x), x=a, x=b, and the x-axis about the x-axis is given by

S A = 2 π ∫ a b f ( x ) 1 + [ f ′ ( x ) ] 2 d x .

(3.23)

Example 3.46

Gabriel's Horn

Gabriel's horn is the solid of revolution obtained by revolving the area of the region bounded by y=1/x and the x-axis for x⩾1 about the x-axis. Show that the surface area of Gabriel's horn is infinite but that its volume is finite.

Solution

After defining f(x)=1/x, we use ParametricPlot3D to visualize a portion of Gabriel's horn in Fig. 3.51.

Figure 3.51 A portion of Gabriel's horn.

f [ x _ ] = 1 / x ;

ParametricPlot3D [ { r , f [ r ] Cos [ t ] , f [ r ] Sin [ t ] } , { r , 1 , 10 } , { t , 0 , 2 Pi } ,

PlotPoints - > { 40 , 40 } , ViewPoint - > { − 1.509 , − 2.739 , 1.294 } ]

Using equation (3.23), the surface area of Gabriel's horn is given by the improper integral

S A = 2 π ∫ 1 ∞ 1 x 1 + 1 x 4 d x = 2 π lim L → ∞ ⁡ ∫ 1 L 1 x 1 + 1 x 4 d x .

f ′ [ x ]

− 1 x 2

step1=Integrate[2Pif[x]Sqrt[1+f′[x]∧2],{x,1,capl}]

ConditionalExpression [ 2 π ( 1 2 ( 2 − ArcSinh [ 1 ] ) + 1 + 1 capl 4 ( − 1 + capl 4 + capl 2 ArcSinh [ capl 2 ] ) 2 1 + capl 4 ) ,

Re [ 1 capl 4 ] ⩾ − 1 & & ( . ( ( ( − 1 + Re [ capl ] ) 4 Im [ capl ] 4 = = − 1 ‖ ( − 1 + Re [ capl ] ) 4 Im [ capl ] 4 ⩾ − 1 ) & & capl ∉ Reals ) ‖ Re [ capl ]

⩾ 0 ) & & ( ( Im [ capl ] + Re [ capl ] = = 0 & & ( 0 < Re [ capl ] < 1 ‖ Im [ capl ] > 0 ‖ Re [ capl ] > 1 ) ) ‖ − 1 < Im [ capl ] < 0 ‖ ( Re [ capl ] ⩾ 0 & & − 1 < Im [ capl ] ⩽ 0 ) ‖ Im [ capl ] < − 1 ) & &

( Im [ capl ] ( Im [ capl ] + Re [ capl ] ) 1 + Im [ capl ] = = 0 ‖ Im [ capl ] 2 ( Im [ capl ] + Re [ capl ] ) 2 ( 1 + Im [ capl ] ) 2 ⩽ 0 ‖ 1 1 + Im [ capl ] = = 1 ‖ 1 1 + Im [ capl ] ⩾ 1 ‖ Im [ capl ] < − 1 ) ]

Limit [ step1 , capl - > Infinity ]

∞

On the other hand, using equation (3.21) the volume of Gabriel's horn is given by the improper integral

V = 2 π ∫ 1 ∞ 1 x 2 d x = π lim L → ∞ ⁡ ∫ 1 L 1 x 2 d x ,

which converges to π.

step1 = Integrate [ Pi f [ x ] ∧ 2 , { x , 1 , capl } ]

ConditionalExpression [ ( 1 − 1 capl ) π , Re [ capl ] > 0 ‖ capl ∉ Reals ]

Limit [ step1 , capl - > Infinity ]

π

Integrate [ Pi f [ x ] ∧ 2 , { x , 1 , Infinity } ]

π  □

3.4 Infinite Sequences and Series

3.4.1 Introduction to Sequences

Sequences and series are usually discussed in the third quarter or second semester of introductory calculus courses. Most students find that it is one of the most difficult topics covered in calculus. A sequence is a function with domain consisting of the positive integers. The terms of the sequence {an} are a1, a2, a3, …. The nth term is an; the (n+1)st term is an+1. If limn→∞⁡an=L, we say that {an} converges to L. If {an} does not converge, {an} diverges. We can sometimes prove that a sequence converges by applying the following theorem.

Theorem 3.7

Every bounded monotonic sequence converges.

A sequence {an} is monotonic if {an} is increasing (an+1⩾an for all n) or decreasing (an+1⩽an for all n).

In particular, Theorem 3.7 gives us the following special cases.

1.  If {an} has positive terms and is eventually decreasing, {an} converges.

2.  If {an} has negative terms and is eventually increasing {an} converges.

After you have defined a sequence, use Table to compute the first few terms of the sequence.

1.  Table[a[n],{n,1,m}] returns the list {a1,a2,a3,…,am}.

2.  Table[a[n],{n,k,m}] returns {ak,ak+1,ak+2,…,am}.

The command ListPlot[listofpoints] plots the list of points listofpoints while ListLinePlot[listofpoints] plots the list of points listofpoints and connects consecutive points with line segments.

Remark 3.7

An extensive database of integer sequences can be found at the On-Line Encyclopedia of Integer Sequences,

http://www.research.att.com/~njas/sequences/Seis.html.

Example 3.47

If an=50nn!, show that limn→∞⁡an=0.

Solution

We remark that the symbol n! in the denominator of an represents the factorial sequence:

n ! = n ⋅ ( n − 1 ) ⋅ ( n − 2 ) ⋅ ⋯ ⋅ 2 ⋅ 1 .

We begin by defining an and then computing the first few terms of the sequence with Table.

Clear [ a ]

a [ n _ ] := 50 ∧ n / n ! ;

afewterms = Table [ a [ n ] , { n , 1 , 10 } ]

{ 50 , 1250 , 62500 3 , 781250 3 , 7812500 3 , 195312500 9 , 9765625000 63 , 61035156250 63 , 3051757812500 567 ,

15258789062500 567 }

N [ afewterms ]

{ 50 . , 1250 . , 20833.3 , 260417 . , 2.60417 × 10 6 , 2.17014 × 10 7 , 1.5501 × 10 8 , 9.68812 × 10 8 ,

5.38229×109,2.69114×1010}

The first few terms increase in magnitude. In fact, this is further confirmed by graphing the first few terms of the sequence with ListPlot in Fig. 3.52 (a). Based on the graph and the values of the first few terms we might incorrectly conclude that the sequence diverges.

Figure 3.52 (a) The first few terms of a_n. (b) The first 75 terms of a_n. (c) The first 75 terms of a_n connected with line segments. (d) Displaying the points together with the first 75 terms of a_n connected with line segments (The Ohio State University colors).

p1 = ListPlot [ afewterms , PlotStyle → { PointSize [ . 025 ] , CMYKColor [ . 03 , 1 , . 63 , . 12 ] } ,

PlotLabel → “(a)” ] ;

However, notice that an+1=50n+1an⇒an+1an=50n+1. Because 50/(n+1)<1 for n>49, we conclude that the sequence is decreasing for n>49. Because it has positive terms, it is bounded below by 0 so the sequence converges by Theorem 3.7. Let L=limn→∞⁡an. Then,

lim n → ∞ ⁡ a n + 1 = lim n → ∞ ⁡ 50 n + 1 a n L = lim n → ∞ ⁡ 50 n + 1 ⋅ L L = 0 .

When we graph a larger number of terms, it is clear that the limit is 0. (See Fig. 3.52 (b).) It is a good exercise to show that for any real value of x, limn→∞⁡xnn!=0.

Use ListLinePlot to connect consecutive points with line segments.

p2 = ListPlot [ Evaluate [ Table [ a [ k ] , { k , 1 , 75 } ] ] ,

PlotStyle → { PointSize [ . 025 ] , CMYKColor [ . 03 , 1 , . 63 , . 12 ] } ,

PlotLabel → “(b)” ] ;

p3 = ListLinePlot [ Evaluate [ Table [ a [ k ] , { k , 1 , 75 } ] ] ,

PlotStyle → CMYKColor [ . 56 , . 47 , . 47 , . 15 ] ,

PlotLabel → “(c)” ] ;

p4 = Show [ p2 , p3 , PlotLabel → “(d)” ] ;

Show[GraphicsGrid[{{p1,p2},{p3,p4}}]] □

Example 3.48

The Rational Numbers and the Calkin–Wilf Sequence

A set S is countably infinite if there is a one-to-one correspondence between the elements of S and the natural numbers. In other words, S is countably infinite if it can be written as a sequence.

Students in introductory set theory classes often learn that the set of rational numbers, Q={nm:n,m≠0 integers} is countably infinite using Cantor's Diagnalization Process.

In 2000, Neil Calkin and Herbert Wilf published a paper introducing the Calkin–Wilf sequence. The Calkin–Wilf sequence is a recursively defined sequence that gives a one-to-one correspondence between the natural numbers and the positive rational numbers. Stating a specific sequence that yields the positive rational numbers proving the countability of the rationals provides an interesting contrast to the traditional method of using Cantor's Diagnalization Process.

The Calkin–Wilf sequence is defined by

q 1 = 1 q i = 1 2 ⌊ q i − 1 ⌋ − q i − 1 + 1  for  i > 1 ,

(3.24)

where ⌊qi−1⌋ represents the integer part of qi−1. The Mathematica command IntegerPart[x] returns the integer part of x. When defining q, observe that we use the form q[i_]:=q[i]= so that Mathematica “remembers” the computed values of qi. Thus, to compute qn, Mathematica need only have computed qn−1 rather than re-computing the previous n−1 qi values. Compute and plot the first n terms of the Calkin–Wilf sequence with n=50 and n=500.

Solution

Clear [ q ]

q [ 1 ] = 1 ;

q [ i _ ] := q [ i ] = 1 / ( 2 IntegerPart [ q [ i − 1 ] ] − q [ i − 1 ] + 1 )

In q1 we compute the first 50 terms of the Calkin–Wilf sequence while in q2 we compute the first 500 terms of the Calkin–Wilf but do not display the result.

q1 = Table [ q [ i ] , { i , 1 , 50 } ]

{ 1 , 1 2 , 2 , 1 3 , 3 2 , 2 3 , 3 , 1 4 , 4 3 , 3 5 , 5 2 , 2 5 , 5 3 , 3 4 , 4 , 1 5 , 5 4 , 4 7 , 7 3 , 3 8 , 8 5 , 5 7 , 7 2 , 2 7 , 7 5 , 5 8 , 8 3 , 3 7 , 7 4 , 4 5 , 5 , 1 6 , 6 5 , 5 9 , 9 4 , 4 11 , 11 7 , 7 10 , 10 3 , 3 11 , 11 8 , 8 13 , 13 5 , 5 12 , 12 7 , 7 9 , 9 2 , 2 9 , 9 7 , 7 12 }

q2 = Table [ q [ i ] , { i , 1 , 500 } ] ;

We then graph both sets of points with ListPlot and show the results in Fig. 3.53.

Figure 3.53 (a) The first 50 terms of the Calkin–Wilf sequence. (b) The first 500 terms of the Calkin–Wilf sequence (Rensselaer Polytechnic Institute colors).

p1 = ListPlot [ q1 , PlotStyle → CMYKColor [ 0 , 1 , 1 , 0 ] ,

PlotLabel → “(a)” ]

p2=ListPlot[q2,PlotStyle→CMYKColor[0,0,0,.61],

PlotLabel → “(b)” ]

Show[GraphicsRow[{p1,p2}]] □

You can use Mathematica to define quite complex sequences. In the following example, we define a sequence that depends on four other sequences. Observe that when defining the terms, we consistently use the form a[n_]:=a[n]=... so that Mathematica “remembers” sequence values computed and need not recompute them to compute subsequent terms.

Example 3.49

The Gauss–Legendre Algorithm

One algorithm to approximate π is the Gauss–Legendre algorithm. Although it is memory intensive, 25 iterations approximate π to around 45 million digits. After defining ak, bk, tk, and pk, π is approximated with π≈(an+bn)24tn. Approximate π to 100 decimal places.

Solution

Remove [ a , b , t , p , pi ]

a [ 0 ] = 1 ;

b [ 0 ] = 1 / Sqrt [ 2 ] ;

t [ 0 ] = 1 / 4 ;

p [ 0 ] = 1 ;

a [ n _ ] := a [ n ] = N [ ( a [ n − 1 ] + b [ n − 1 ] ) / 2 , 10000 ] ;

b [ n _ ] := b [ n ] = N [ Sqrt [ a [ n − 1 ] b [ n − 1 ] ] , 10000 ] ;

t [ n _ ] := t [ n ] = N [ t [ n − 1 ] − p [ n − 1 ] ( a [ n − 1 ] − a [ n ] ) ∧ 2 , 10000 ] ;

p [ n _ ] := p [ n ] = N [ 2 p [ n − 1 ] , 10000 ] ;

π is then approximated with

pi[n_]:=N[(a[n]+b[n])∧2/(4t[n]),100]

To increase the accuracy of your approximation, increase 10000 and 100 to your desired values. However, even on relatively fast machines, if you replace all these numbers with a number such as 1000000, the computation time will increase considerably.

For illustrate purposes, we approximate π with π≈(a6+b6)24t6. Observe that the result is accurate to at least 100 decimal places.

pi [ 6 ]

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

N [ Pi , 100 ]

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 □

3.4.2 Introduction to Infinite Series

An infinite series is a series of the form

∑ k = 1 ∞ a k = a 1 + a 2 + a 3 + ⋯

(3.25)

where {an} is a sequence. The nth partial sum of (3.25) is

s n = ∑ k = 1 n a k = a 1 + a 2 + ⋯ + a n .

(3.26)

Notice that the partial sums of the series (3.25) form a sequence {sn}. Hence, we say that the infinite series (3.25) converges to L if the sequence of partial sums {sn} converges to L and write

∑ k = 1 ∞ a k = L .

The infinite series (3.25) diverges if the sequence of partial sums diverges. Given the infinite series (3.25),

Sum[a[k],{k,1,n}]

calculates the nth partial sum (3.26). In some cases, if the infinite series (3.25) converges,

Sum[a[k],{k,1,Infinity}]

can compute the value of the infinite sum. In addition to using Sum to compute finite and infinite sums, you can use the button on the Basic Math Input palette to calculate sums. You should think of the Sum function as a “fragile” command and be certain to carefully examine its results and double check their validity.

Example 3.50

Determine whether each series converges or diverges. If the series converges, find its sum. (a) ∑k=1∞(−1)k+1; (b) ∑k=2∞2k2−1; and (c) ∑k=0∞ark.

Solution

For (a), we compute the nth partial sum (3.26) in sn with Sum.

sn = Sum [ ( − 1 ) ∧ ( k + 1 ) , { k , 1 , n } ]

1 2 ( 1 + ( − 1 ) 1 + n )

Notice that the odd partial sums are 1: s2n+1=12((−1)2n+1+1+1)=12(1+1)=1 while the even partial sums are 0: s2n=12((−1)2n+1+1)=12(−1+1)=0. We confirm that the limit of the partial sums does not exist with Limit. Mathematica's result indicates that it cannot determine the limit. The series diverges.

Limit [ sn , n → Infinity ]

1 2 ( 1 + e 2 i Interval [ { 0 , π } ] )

Similarly, when we attempt to compute the infinite sum with Sum, Mathematica is able to determine that the partial sums diverge, which means that the infinite series diverges.

For (b), we have a telescoping series. Using partial fractions,

∑ k = 2 ∞ 2 k 2 − 1 = ∑ k = 2 ∞ ( 1 k − 1 − 1 k + 1 ) = ( 1 − 1 3 ) + ( 1 2 − 1 4 ) + ( 1 3 − 1 5 ) + ⋯ + ( 1 n − 2 − 1 n ) + ( 1 n − 1 − 1 n + 1 ) + …

we see that the nth partial sum is given by

s n = 3 2 − 1 n − 1 n + 1

and sn→3/2 as n→∞ so the series converges to 3/2:

∑ k = 2 ∞ 2 k 2 − 1 = 3 2 .

We perform the same steps with Mathematica using Sum, Apart, and Limit.

Apart computes the partial fraction decomposition of a rational expression.

sn = Sum [ 1 / ( k − 1 ) − 1 / ( k + 1 ) , { k , 2 , n } ]

−2−n+3n22n(1+n)

Apart [ sn ]

3 2 − 1 n − 1 1 + n

Limit [ sn , n → Infinity ]

3 2

(c) A series of the form ∑k=0∞ark is called a geometric series. We compute the nth partial sum of the geometric series with Sum.

sn = Sum [ a r ∧ k , { k , 0 , n } ]

a ( − 1 + r 1 + n ) − 1 + r

When using Limit to determine the limit of sn as n→∞, we see that Mathematica returns the limit unevaluated because Mathematica does not know the value of r.

Limit [ sn , n → Infinity ]

Limit [ a ( − 1 + r 1 + n ) − 1 + r , n → ∞ ]

In fact, the geometric series diverges if |r|⩾1 and converges if |r|<1. Observe that if we simply compute the sum with Sum, Mathematica returns a/(1−r) which is correct if |r|<1 but incorrect if |r|⩾1.

Sum [ a r ∧ k , { k , 0 , Infinity } ]

a 1 − r

However, the result of entering

is correct because the series ∑k=0∞(−53)k is geometric with |r|=5/3⩾1 and, consequently, diverges. Similarly,

Sum [ 9 ( 1 / 10 ) ∧ k , { k , 1 , Infinity } ]

1

is correct because ∑k=1∞9(110)k is geometric with a=9/10 and r=1/10 so the series 0.999… converges to

a 1 − r = 9 / 10 1 − 1 / 10 = 1 .

 □

3.4.3 Convergence Tests

Frequently used convergence tests for infinite series are stated in the following theorems. Note that the infinite series ∑k=1∞ak converges absolutely means that ∑k=1∞|ak| converges. If an infinite series converges absolutely, it converges. If a series converges, but does not converge absolutely, we say that the series conditionally converges.

The alternating harmonic series, ∑k=1∞(−1)k+11/k is an example of a series that converges conditionally. It is a good exercise to show that the alternating harmonic series converges to ln⁡2. On the other hand, the harmonic series, ∑k=1∞1/k diverges to +∞.

Theorem 3.8

The Divergence Test

Let ∑k=1∞ak be an infinite series. If limk→∞⁡ak≠0, then ∑k=1∞ak diverges.

Theorem 3.9

The Integral Test

Let ∑k=1∞ak be an infinite series with positive terms. If f(x) is a decreasing continuous function for which f(k)=ak for all k, then ∑k=1∞ak and ∫1∞f(x)dx either both converge or both diverge.

Theorem 3.10

The Ratio Test

Let ∑k=1∞ak be an infinite series and let ρ=limk→∞⁡|ak+1ak|.

1.  If ρ<1, ∑k=1∞ak converges absolutely.

2.  If ρ>1, ∑k=1∞ak diverges.

3.  If ρ=1, the Ratio Test is inconclusive.

Theorem 3.11

The Root Test

Let ∑k=1∞ak be an infinite series and let ρ=limk→∞⁡|ak|k.

1.  If ρ<1, ∑k=1∞ak converges absolutely.

2.  If ρ>1, ∑k=1∞ak diverges.

3.  If ρ=1, the Root Test is inconclusive.

Theorem 3.12

The Limit Comparison Test

Let ∑k=1∞ak and ∑k=1∞bk be infinite series with positive terms and let L=limk→∞⁡akbk. If 0<L<∞, then either both series converge or both series diverge.

Example 3.51

Determine whether each series converges or diverges. (a) ∑k=1∞(1+1k)k; (b) ∑k=1∞1kp; (c) ∑k=1∞k3k; (d) ∑k=1∞(k!)2(2k)!; (e) ∑k=1∞(k4k+1)k; (f) ∑k=1∞2k+1(k+1)(2k+1).

Solution

(a) Using Limit, we see that the limit of the terms is e≠0 so the series diverges by the Divergence Test, Theorem 3.8.

Limit [ ( 1 + 1 / k ) ∧ k , k → Infinity ]

e

It is a very good exercise to show that the limit of the terms of the series is e by hand. Let L=limk→∞⁡(1+1/k)k. Take the logarithm of each side of this equation and apply L'Hôpital's rule:

ln ⁡ L = lim k → ∞ ⁡ ln ⁡ ( 1 + 1 k ) k ln ⁡ L = lim k → ∞ ⁡ k ln ⁡ ( 1 + 1 k ) ln ⁡ L = lim k → ∞ ⁡ ln ⁡ ( 1 + 1 k ) 1 k ln ⁡ L = lim k → ∞ ⁡ 1 1 + 1 k ⋅ − 1 k 2 − 1 k 2 ln ⁡ L = 1 .

Exponentiating yields L=eln⁡L=e1=e.

(b) A series of the form ∑k=1∞1kp (p>0) is called a p-series. Let f(x)=x−p. Then, f(x) is continuous and decreasing for x⩾1, f(k)=k−p and

∫ 1 ∞ x − p d x = { ∞ , if  p ⩽ 1 1 / ( p − 1 ) , if  p > 1

so the p-series converges if p>1 and diverges if p⩽1. If p=1, the series ∑k=1∞1k is called the harmonic series. Observe that Mathematica's result is given in terms of a “conditional expression.” In this example, the interpretation is that ∫1∞x−pdx=1/(p−1) provided that the real part of p is greater than 1.

Clear [ x , p ]

s1 = Integrate [ x ∧ ( − p ) , { x , 1 , Infinity } ]

ConditionalExpression [ 1 − 1 + p , Re [ p ] > 1 ]

(c) Let f(x)=x⋅3−x. Then, f(k)=k⋅3−k and f(x) is decreasing for x>1/ln⁡3 because f′(x)<0 for x>1/ln⁡3.

f [ x _ ] = x 3 ∧ ( − x ) ;

Factor [ f ′ [ x ] ]

− 3 − x ( − 1 + x Log [ 3 ] )

Solve [ − 1 + x Log [ 3 ] = = 0 , x ]

{{x→1Log[3]}}

Using Integrate, we see that the improper integral ∫1∞f(x)dx converges.

ival = Integrate [ f [ x ] , { x , 1 , Infinity } ]

N [ ival ]

1 + Log [ 3 ] 3 Log [ 3 ] 2

0.579592

Thus, by the Integral Test, Theorem 3.9, we conclude that the series converges. Note that when applying the Integral Test, if the improper integral converges its value is not the value of the sum of the series. In this case, we see that Mathematica is able to evaluate the sum with Sum and the series converges to 3/4.

Sum [ k 3 ∧ ( − k ) , { k , 1 , Infinity } ]

3 4

(d) If ak contains factorial functions, the Ratio Test is often a good first test to try. After defining ak we compute

lim k → ∞ ⁡ a k + 1 a k = lim k → ∞ ⁡ [ ( k + 1 ) ! ] 2 [ 2 ( k + 1 ) ] ( k ! ) 2 ( 2 k ) ! = lim k → ∞ ⁡ ( k + 1 ) ! ⋅ ( k + 1 ) ! k ! ⋅ k ! ( 2 k ) ! ( 2 k + 2 ) ! = lim k → ∞ ⁡ ( k + 1 ) 2 ( 2 k + 2 ) ( 2 k + 1 ) = lim k → ∞ ⁡ ( k + 1 ) 2 ( 2 k + 1 ) = 1 4 .

Because 1/4<1, the series converges by the Ratio Test. We confirm these results with Mathematica.

Remark 3.8

Use FullSimplify instead of Simplify to simplify expressions involving factorials.

a [ k _ ] = ( k ! ) ∧ 2 / ( 2 k ) ! ;

s1 = FullSimplify [ a [ k + 1 ] / a [ k ] ]

1 + k 2 + 4 k

Limit [ s1 , k → Infinity ]

1 4

We illustrate that we can evaluate the sum using Sum and approximate it with N as follows.

ev=Sum[a[k],{k,1,Infinity}]

1 27 ( 9 + 2 3 π )

N [ ev ]

0.7364

(e) Because

lim k → ∞ ⁡ ( k 4 k + 1 ) k k = lim k → ∞ ⁡ k 4 k + 1 = 1 4 < 1 ,

the series converges by the Root Test.

a [ k _ ] = ( k / ( 4 k + 1 ) ) ∧ k ;

Limit [ a [ k ] ∧ ( 1 / k ) , k → Infinity ]

1 4

As with (d), we can approximate the sum with N and Sum.

ev = Sum [ a [ k ] , { k , 1 , Infinity } ]

∑ k = 1 ∞ ( k 1 + 4 k ) k

N [ ev ]

0.265757

(f) We use the Limit Comparison Test and compare the series to ∑k=1∞kkk=∑k=1∞1k, which diverges because it is a p-series with p=1. Because

0 < lim k → ∞ ⁡ 2 k + 1 ( k + 1 ) ( 2 k + 1 ) 1 k = 1 < ∞

and the harmonic series diverges, the series diverges by the Limit Comparison Test.

a [ k _ ] = ( 2 Sqrt [ k ] + 1 ) / ( ( Sqrt [ k ] + 1 ) ( 2 k + 1 ) ) ;

b [ k _ ] = 1 / k ;

Limit [ a [ k ] / b [ k ] , k → Infinity ]

1  □

3.4.4 Alternating Series

An alternating series is a series of the form

∑ k = 1 ∞ ( − 1 ) k a k or ∑ k = 1 ∞ ( − 1 ) k + 1 a k

(3.27)

where {ak} is a sequence with positive terms.

Theorem 3.13

Alternating Series Test

If {ak} is decreasing and limk→∞⁡ak=0, the alternating series (3.27) converges.

Definition 3.3

The alternating series (3.27) converges absolutely if ∑k=1∞ak converges.

Theorem 3.14

If the alternating series (3.27) converges absolutely, it converges.

Definition 3.4

If the alternating series (3.27) converges but does not converge absolutely, we say that the alternating series, (3.27), conditionally converges.

Example 3.52

Determine whether each series converges or diverges. If the series converges, determine whether the convergence is conditional or absolute. (a) ∑k=1∞(−1)k+1k; (b) ∑k=1∞(−1)k+1(k+1)!4k(k!)2; (c) ∑k=1∞(−1)k+1(1+1k)k.

Solution

(a) Because {1/k} is decreasing and 1/k→0 as k→∞, the series converges. The series does not converge absolutely because the harmonic series diverges. Hence, ∑k=1∞(−1)k+1k, which is called the alternating harmonic series, converges conditionally. We see that this series converges to ln⁡2 with Sum.

a [ k _ ] = ( − 1 ) ∧ ( k + 1 ) / k ;

Sum [ a [ k ] , { k , 1 , Infinity } ]

Log [ 2 ]

(b) We test for absolute convergence first using the Ratio Test. Because

lim k → ∞ ⁡ ( ( k + 1 ) + 1 ) ! 4 k + 1 [ ( k + 1 ) ! ] 2 ( k + 1 ) ! 4 k ( k ! ) 2 = lim k → ∞ ⁡ k + 2 4 ( k + 1 ) 2 = 0 < 1 ,

a [ k _ ] = ( k + 1 ) ! / ( 4 ∧ k ( k ! ) ∧ 2 ) ;

s1 = FullSimplify [ a [ k + 1 ] / a [ k ] ]

Limit [ s1 , k → Infinity ]

2 + k 4 ( 1 + k ) 2

0

the series converges absolutely by the Ratio Test. Absolute convergence implies convergence so the series converges.

(c) Because limk→∞⁡(1+1k)k=e, limk→∞⁡(−1)k+1(1+1k)k does not exist, so the series diverges by the Divergence Test. We confirm that the limit of the terms is not zero with Limit.

Sum [ ( − 1 ) ∧ ( k + 1 ) a [ k ] , { k , 1 , Infinity } ]

− 3 − 4 e 1 / 4 4 e 1 / 4

Limit [ a [ k ] , k → Infinity ]

e1+2iInterval[{0,π}] □

3.4.5 Power Series

Let x0 be a number. A power series in x−x0 is a series of the form

∑ k = 0 ∞ a k ( x − x 0 ) k .

(3.28)

A fundamental problem is determining the values of x, if any, for which the power series converges, the interval of convergence.

Theorem 3.15

For the power series (3.28), exactly one of the following is true.

1.  The power series converges absolutely for all values of x. The interval of convergence is (−∞,∞).

2.  There is a positive number r so that the series converges absolutely if x0−r<x<x0+r. The series may or may not converge at x=x0−r and x=x0+r. The interval of convergence will be one of (x0−r,x0+r), [x0−r,x0+r), (x0−r,x0+r], or [x0−r,x0+r].

3.  The series converges only if x=x0. The interval of convergence is {x0}.

Example 3.53

Determine the interval of convergence for each of the following power series. (a) ∑k=0∞(−1)k(2k+1)!x2k+1; (b) ∑k=0∞k!1000k(x−1)k; (c) ∑k=1∞2kk(x−4)k.

Solution

(a) We test for absolute convergence first using the Ratio Test. Because

lim k → ∞ ⁡ | ( − 1 ) k + 1 ( 2 ( k + 1 ) + 1 ) ! x 2 ( k + 1 ) + 1 ( − 1 ) k ( 2 k + 1 ) ! x 2 k + 1 | = lim k → ∞ ⁡ 1 2 ( k + 1 ) ( 2 k + 3 ) x 2 = 0 < 1

a [ x _ , k _ ] = ( − 1 ) ∧ k / ( 2 k + 1 ) ! x ∧ ( 2 k + 1 ) ;

s1 = FullSimplify [ a [ x , k + 1 ] / a [ x , k ] ]

Limit [ s1 , k → Infinity ]

− x 2 6 + 10 k + 4 k 2

0

for all values of x, we conclude that the series converges absolutely for all values of x; the interval of convergence is (−∞,∞). In fact, we will see later that this series converges to sin⁡x:

sin ⁡ x = ∑ k = 0 ∞ ( − 1 ) k + 1 ( 2 k + 1 ) ! x 2 k + 1 = x − 1 3 ! x 3 + 1 5 ! x 5 − 1 7 ! x 7 + … ,

which means that the partial sums of the series converge to sin⁡x. Graphically, we can visualize this by graphing partial sums of the series together with the graph of y=sin⁡x. Note that the partial sums of a series are a recursively defined function: sn=sn−1+an, s0=a0. We use this observation to define p to be the nth partial sum of the series. We use the form p[x_,n_]:=p[x,n]=... so that Mathematica “remembers” the partial sums computed. That is, once p[x,3] is computed, Mathematica need not recompute p[x,3] when computing p[x,4].

In Fig. 3.54 we graph pn(x)=∑k=0n(−1)k(2k+1)!x2k+1 together with y=sin⁡x for n=1, 5, and 10. In the graphs, notice that as n increases, the graphs of pn(x) more closely resemble the graph of y=sin⁡x.

Figure 3.54 y=sin⁡x together with the graphs of p₁(x), p₅(x), and p₁₀(x).

When you use Tooltip, placing the cursor over the plot shows you the function being plotted.

Clear [ p ]

p [ x _ , 0 ] = a [ x , 0 ] ;

p[x_,n_]:=p[x,n]=p[x,n−1]+a[x,n]

p [ x , 2 ]

x − x 3 6 + x 5 120

p1 = Plot [ Tooltip [ { Sin [ x ] , p [ x , 1 ] , p [ x , 5 ] } ] , { x , − 2 Pi , 2 Pi } ,

PlotStyle → { { CMYKColor [ 1 , . 7 , . 07 , . 45 ] } , { CMYKColor [ 0 , . 63 , 1 , 0 ] } ,

{ CMYKColor [ . 96 , . 09 , . 32 , . 29 ] } } ,

PlotRange → { − Pi , Pi } , AspectRatio → Automatic ]

We use Manipulate to investigate how n affects the situation with (see Fig. 3.55)

Figure 3.55 Using Manipulate to investigate the situation.

p2 = Manipulate [ Plot [ Tooltip [ { Sin [ x ] , p [ x , n ] } ] ,

{ x , − 2 Pi , 2 Pi } ,

PlotStyle → { { CMYKColor [ 1 , . 7 , . 07 , . 45 ] } , { CMYKColor [ 0 , . 63 , 1 , 0 ] } } ,

PlotRange → { − 2 Pi , 2 Pi } , AspectRatio → Automatic ] , { { n , 5 } , 1 , 25 , 1 } ]

(b) As in (a), we test for absolute convergence first using the Ratio Test:

lim k → ∞ ⁡ | ( k + 1 ) k ! 1000 k + 1 ( x − 1 ) k + 1 k ! 1000 k ( x − 1 ) k | = 1 1000 ( k + 1 ) | x − 1 | = { 0 , if  x = 1 ∞ , if  x ≠ 1 .

a [ x _ , k _ ] = k ! / 1000 ∧ k ( x − 1 ) ∧ k ;

s1 = FullSimplify [ a [ x , k + 1 ] / a [ x , k ] ]

Limit[s1,k→Infinity]

( 1 + k ) ( − 1 + x ) 1000

( − 1 + x ) ∞

Be careful of your interpretation of the result of the Limit command because Mathematica does not consider the case x=1 separately: if x=1 the limit is 0. Because 0<1 the series converges by the Ratio Test.

The series converges only if x=1; the interval of convergence is {1}. You should observe that if you graph several partial sums for “small” values of n, you might incorrectly conclude that the series converges. (c) Use the Ratio Test to check absolute convergence first:

lim k → ∞ ⁡ | 2 k + 1 k + 1 ( x − 4 ) k + 1 2 k k ( x − 4 ) k | = lim k → ∞ ⁡ 2 k k + 1 | x − 4 | = 2 | x − 4 | .

By the Ratio Test, the series converges absolutely if 2|x−4|<1. We solve this inequality for x with Reduce to see that 2|x−4|<1 if 7/2<x<9/2.

Clear [ a , s1 , k ]

a [ x _ , k _ ] = 2 ∧ k / Sqrt [ k ] ( x − 4 ) ∧ k ;

s1 = Simplify [ Abs [ a [ x , k + 1 ] / a [ x , k ] ] ]

Limit [ s1 , k → Infinity ]

2 Abs [ k 1 + k ( − 4 + x ) ]

2 Abs [ − 4 + x ]

Use Reduce to solve the inequality.

Reduce [ 2 Abs [ x − 4 ] < 1 , x ]

7 2 < Re [ x ] < 9 2 & & − 1 2 − 63 + 32 Re [ x ] − 4 Re [ x ] 2 < Im [ x ]

< 1 2 − 63 + 32 Re [ x ] − 4 Re [ x ] 2

From the output, we see that for real values of x, the inequality is satisfied for 7/2<x<9/2. We check x=7/2 and x=9/2 separately. If x=7/2, the series becomes ∑k=1∞(−1)k1k, which converges conditionally.

Simplify [ a [ x , k ] /. x → 7 / 2 ]

( − 1 ) k k

On the other hand, if x=9/2,

Simplify [ a [ x , k ] /. x → 9 / 2 ]

1 k

the series is ∑k=1∞1k, which diverges. We conclude that the interval of convergence is [7/2,9/2). □

3.4.6 Taylor and Maclaurin Series

Let y=f(x) be a function with derivatives of all orders at x=x0. The Taylor series for f(x) about x=x0 is

∑ k = 0 ∞ f ( k ) ( x 0 ) k ! ( x − x 0 ) k .

(3.29)

The Maclaurin series for f(x) is the Taylor series for f(x) about x=0. If y=f(x) has derivatives up to at least order n at x=x0, the nth degree Taylor polynomial for f(x) about x=x0 is

p n ( x ) = ∑ k = 0 n f ( k ) ( x 0 ) k ! ( x − x 0 ) k .

(3.30)

The nth degree Maclaurin polynomial for f(x) is the nth degree Taylor polynomial for f(x) about x=0. Generally, finding Taylor and Maclaurin series using the definition is a tedious task at best.

Example 3.54

Find the first few terms of (a) the Maclaurin series and (b) the Taylor series about x=π/4 for f(x)=tan⁡x.

Solution

(a) After defining f(x)=tan⁡x, we use Table together with /. and D to compute f(k)(0)/k! for k=0, 1, …, 8.

Use Short to obtain an abbreviated result. Many terms will be missing, but with Short, you will see the beginning and end of your result.

Clear [ f ]

f [ x _ ] = Tan [ x ] ;

t1 = Table [ { k , D [ f [ x ] , { x , k } ] , D [ f [ x ] , { x , k } ] /. x → 0 } , { k , 0 , 8 } ] ;

Short [ t1 ]

{ { 0 , Tan [ x ] , 0 } , 〈 〈 7 〉 〉 , { 8 , 7936 Sec [ x ] 8 Tan [ x ] + 24576 〈 〈 1 〉 〉 6 〈 〈 1 〉 〉 3 + 〈 〈 1 〉 〉 + 128 Sec [ x ] 2

Tan [ x ] 7 , 0 } }

To see these results in tabular form, enter

t1 // TableForm

For length considerations, the resulting output is not shown here. Another way of approaching the problem is to use Manipulate. See Fig. 3.56.

Figure 3.56 With Manipulate, we can adjust the function and function values.

Manipulate [ { k , D [ f [ x ] , { x , k } ] // FullSimplify , D [ f [ x ] , { x , k } ] /. x → 0 ,

D [ f [ x ] , { x , k } ] /. x → 0 // N } , { { k , 5 } , 0 , 25 , 1 } ]

Using the values in the table or from the Manipulate object, we apply the definition to see that the Maclaurin series is

∑ k = 0 ∞ f ( k ) ( 0 ) k ! x k = x + 1 3 x 3 + 2 15 x 5 + 17 315 x 7 + …

For (b), we repeat (a) using x=π/4 instead of x=0

Manipulate [ { k , D [ f [ x ] , { x , k } ] // FullSimplify , D [ f [ x ] , { x , k } ] /. x → Pi / 4 ,

D [ f [ x ] , { x , k } ] /. x → Pi / 4 // N } , { { k , 1 } , 0 , 25 , 1 } ]

and then apply the definition to see that the Taylor series about x=π/4 is

∑ k = 0 ∞ f ( k ) ( x 0 ) k ! ( x − x 0 ) k = 1 + 2 ( x − π 4 ) + 2 ( x − π 4 ) 2 + 8 3 ( x − π 4 ) 3 + 10 3 ( x − π 4 ) 4 + 64 15 ( x − π 4 ) 5 + 244 45 ( x − π 4 ) 6 + …

From the series, we can see various Taylor and Maclaurin polynomials. For example, the third Maclaurin polynomial is

p 3 ( x ) = x + 1 3 x 3

and the 4th degree Taylor polynomial about x=π/4 is

p 4 ( x ) = 1 + 2 ( x − π 4 ) + 2 ( x − π 4 ) 2 + 8 3 ( x − π 4 ) 3 + 10 3 ( x − π 4 ) 4 .

See Fig. 3.57. □

Figure 3.57 We use Manipulate to investigate series for the tangent function.

The command Series[f[x],{x,x0,n}] computes (3.29) to (at least) order n−1. Because of the O-term in the result that represents the terms that are omitted from the power series for f(x) expanded about the point x=x0, the result of entering a Series command is not a function that can be evaluated if x is a particular number. We remove the remainder (O-) term of the power series Series[f[x],{x,x0,n}] with the command Normal and can then evaluate the resulting polynomial for particular values of x.

Example 3.55

Find the first few terms of the Taylor series for f(x) about x=x0. (a) f(x)=cos⁡x, x=0; (b) f(x)=1/x2, x=1.

Solution

Entering

Series [ Cos [ x ] , { x , 0 , 4 } ]

1 − x 2 2 + x 4 24 + O [ x ] 5

computes the Maclaurin series to order 4. Entering

Series [ Cos [ x ] , { x , 0 , 14 } ]

1 − x 2 2 + x 4 24 − x 6 720 + x 8 40320 − x 10 3628800 + x 12 479001600 − x 14 87178291200 + O [ x ] 15

computes the Maclaurin series to order 14. In this case, the Maclaurin series for cos⁡x converges to cos⁡x for all real x. To graphically see this, we define the function p. Given n, p[n] returns the Maclaurin polynomial of degree n for cos⁡x.

p [ n _ ] := Series [ Cos [ x ] , { x , 0 , n } ] // Normal

p [ 3 ]

1 − x 2 2

We then graph cos⁡x together with the Maclaurin polynomial of degree n=2, 4, 8, and 16 on the interval [−3π/2,3π/2] in Fig. 3.58. Notice that as n increases, the graph of the Maclaurin polynomial more closely resembles the graph of cos⁡x. We would see the same pattern if we increased the length of the interval and the value of n.

Figure 3.58 Using Manipulate to investigate graphs of y=cos⁡x together with plots of several of its Maclaurin polynomials (University of Pennsylvania colors).

Manipulate [ Plot [ Evaluate [ Tooltip [ { Cos [ x ] , p [ n ] } ] ] , { x , − 3 Pi / 2 , 3 Pi / 2 } ,

PlotStyle → { { CMYKColor [ 0 , 1 , . 65 , . 34 ] } , { CMYKColor [ 1 , . 65 , 0 , . 3 ] } } ,

PlotRange → { { − 3 Pi / 2 , 3 Pi / 2 } , { − 3 Pi / 2 , 3 Pi / 2 } } , AspectRatio → Automatic ] ,

{ { n , 4 } , { 2 , 4 , 8 , 16 , 32 , 64 } , ControlType → Setter } ]

(b) After defining f(x)=1/x2, we compute the first 10 terms of the Taylor series for f(x) about x=1 with Series.

f [ x _ ] = 1 / x ∧ 2 ;

p10 = Series [ f [ x ] , { x , 1 , 10 } ]

1−2(x−1)+3(x−1)2−4(x−1)3+5(x−1)4−6(x−1)5+7(x−1)6−8(x−1)7+9(x−1)8−10(x−1)9+11(x−1)10+O[x−1]11

In this case, the pattern for the series is relatively easy to see: the Taylor series for f(x) about x=1 is

∑ k = 0 ∞ ( − 1 ) k ( k + 1 ) ( x − 1 ) k .

This series converges absolutely if

lim k → ∞ ⁡ | ( − 1 ) k + 1 ( k + 2 ) ( x − 1 ) k + 1 ( − 1 ) k ( k + 1 ) ( x − 1 ) k | = | x − 1 | < 1

or 0<x<2. The series diverges if x=0 and x=2. In this case, the series converges to f(x) on the interval (0,2).

a [ x _ , k _ ] = ( − 1 ) ∧ k ( k + 1 ) ( x − 1 ) ∧ k ;

s1 = FullSimplify [ Abs [ a [ x , k + 1 ] / a [ x , k ] ] ]

Abs [ ( 2 + k ) ( − 1 + x ) 1 + k ]

s2 = Limit [ s1 , k → Infinity ]

Abs [ − 1 + x ]

Reduce [ s2 < 1 , x ]

0 < Re [ x ] < 2 & & − 2 Re [ x ] − Re [ x ] 2 < Im [ x ] < 2 Re [ x ] − Re [ x ] 2

To see this, we use Manipulate to graph f(x) together with the Taylor polynomial for f(x) about x=1 of degree n for large n. Regardless of the size of n, the graphs of f(x) and the Taylor polynomial closely resemble each other on the interval (0,2) – but not at the endpoints or outside the interval. (See Fig. 3.59.)

Figure 3.59 Graphs of f(x) together with the various Taylor polynomials about x = 1.

p [ n _ ] := Series [ f [ x ] , { x , 1 , n } ] // Normal

Manipulate [ Plot [ Evaluate [ Tooltip [ { f [ x ] , p [ n ] } ] ] ,

{ x , 0 , 3 } , PlotStyle → { { CMYKColor [ 0 , 1 , . 65 , . 34 ] } , { CMYKColor [ 1 , . 65 , 0 , . 3 ] } } ,

PlotRange → { { 0 , 3 } , { − 3 / 2 , 3 / 2 } } , AspectRatio → Automatic ] ,

{{n,10},{5,10,15,25,30,35,40,45,50,55,60},ControlType→Setter}] □

3.4.7 Taylor's Theorem

Taylor's theorem states the relationship between f(x) and the Taylor series for f(x) about x=x0.

Theorem 3.16

Taylor's Theorem

Let y=f(x) have (at least) n+1 derivatives on an interval I containing x=x0. Then, for every number x∈I, there is a number z between x and x0 so that

f ( x ) = p n ( x ) + R n ( x ) ,

where pn(x) is given by equation (3.30) and

R n ( x ) = f ( n + 1 ) ( z ) ( n + 1 ) ! ( x − x 0 ) n + 1 .

(3.31)

Example 3.56

Use Taylor's theorem to show that

sin ⁡ x = ∑ k = 0 ∞ ( − 1 ) k ( 2 k + 1 ) ! x 2 k + 1

Solution

Let f(x)=sin⁡x. Then, for each value of x, there is a number z between 0 and x so that sin⁡x=pn(x)+Rn(x) where pn(x)=∑k=0nf(k)(0)k!xk and Rn(x)=f(n+1)(z)(n+1)!xn+1. Regardless of the value of n, f(n+1)(z) is one of sin⁡z, −sin⁡z, cos⁡z, or −cos⁡z, which are all bounded by 1. Then,

| sin ⁡ x − p n ( x ) | = | f ( n + 1 ) ( z ) ( n + 1 ) ! x n + 1 | | sin ⁡ x − p n ( x ) | ⩽ 1 ( n + 1 ) ! | x | n + 1

and xn/n!→0 as n→∞ for all real values of x.

You should remember that the number z in Rn(x) is guaranteed to exist by Taylor's theorem. However, from a practical point of view, you would rarely (if ever) need to compute the z value for a particular x value.

For illustrative purposes, we show the difficulties. Suppose we wish to approximate sin⁡π/180 using the Maclaurin polynomial of degree 4, p4(x)=x−16x3, for sin⁡x. The fourth remainder is R4(x)=1120cos⁡zx5.

The Maclaurin polynomial of degree 4 for sin⁡x is ∑k=04f(k)(0)k!x4= 0+x+0⋅x2+−13!x3+0⋅x4.

Clear [ f ]

f [ x _ ] = Sin [ x ] ;

r5 = D [ f [ z ] , { z , 5 } ] / 5 ! x ∧ 5

1 120 x 5 Cos [ z ]

If x=π/180 there is a number z between 0 and π/180 so that

| R 4 ( π 180 ) | = 1 120 cos ⁡ z ( π 180 ) 5 ⩽ 1 120 ( π 180 ) 5 ≈ 0.135 × 10 − 10 ,

which shows us that the maximum the error can be is 1120(π180)5≈0.135×10−10.

maxerror = N [ 1 / 120 ⁎ ( Pi / 180 ) ∧ 5 ]

1.349601623163255 ` * − ∧ 11

Abstractly, the exact error can be computed. By Taylor's theorem, z satisfies

f ( π 180 ) = p 4 ( π 180 ) + R 4 ( π 180 ) sin ⁡ π 180 = 1 180 π − 1 34992000 π 3 + 1 22674816000000 π 5 cos ⁡ z 0 = 1 180 π − 1 34992000 π 3 + 1 22674816000000 π 5 cos ⁡ z − sin ⁡ π 180 .

We graph the right-hand side of this equation with Plot in Fig. 3.60. The exact value of z is the z-coordinate of the point where the graph intersects the z-axis.

Figure 3.60 Finding z (Observe that Mathematica's default color is similar to Columbia University's blue).

p4 = Series [ f [ x ] , { x , 0 , 4 } ] // Normal

x − x 3 6

exval = Sin [ Pi / 180 ]

p4b = p4 /. x → Pi / 180

r5b = r5 /. x → Pi / 180

Sin [ π 180 ]

π 180 − π 3 34992000

π 5 Cos [ z ] 22674816000000

toplot = r5b + p4b − exval ;

Plot [ toplot , { z , 0 , Pi / 180 } ]

We can use FindRoot to approximate z, if we increase the number of digits carried in floating point calculations with WorkingPrecision.

exz = FindRoot [ toplot = = 0 , { z , 0 , . 004 } , WorkingPrecision → 32 ]

{ z → 0.0038086149165541606417429516417308 }

Alternatively, we can compute the exact value of z with Solve

cz = Solve [ toplot = = 0 , z ]

{ { z → ConditionalExpression [ − ArcCos [ 648000 ( − 194400 π + π 3 + 34992000 Sin [ π 180 ] ) π 5 ] + 2 π C [ 1 ] ,

C[1]∈Integers]},{z→ConditionalExpression[ArcCos[648000(−194400π+π3+34992000Sin[π180])π5]

+ 2 π C [ 1 ] , C [ 1 ] ∈ Integers ] } }

and then approximate the result with N.

N [ cz ]

{ { z → ConditionalExpression [ − 0.00384232 + 6.28319 C [ 1 ] , C [ 1 ] ∈ Integers ] } ,

{z→ConditionalExpression[0.00384232+6.28319C[1],C[1]∈Integers]}} □

Example 3.57

Newton's Approximation of π.

Newton's Approximation of π is a fascinating combination of calculus, infinite series, and geometry. Newton began the approximation after proving the Binomial Theorem.

Theorem 3.17

Binomial Theorem

For −1<x<1,

( 1 + x ) r = 1 + ∑ k = 1 ∞ ( r k ) x k ,

(3.32)

where (r1)=1, (r2)=r(r−1)2!, and (rk)=r(r−1)(r−2)⋯(r−k+1)k! for k⩾3.

We use Series to find the first few terms of the Binomial series.

Series [ ( 1 + x ) ∧ r , { x , 0 , 5 } ]

1 + r x + 1 2 ( − 1 + r ) r x 2 + 1 6 ( − 2 + r ) ( − 1 + r ) r x 3 + 1 24 ( − 3 + r ) ( − 2 + r ) ( − 1 + r ) r x 4 + 1 120 ( − 4 + r ) ( − 3 + r ) ( − 2 + r ) ( − 1 + r ) r x 5 + O [ x ] 6

Newton then used the Binomial series to find the first few terms of the series for 1−x=(1−x)1/2.

f [ x _ ] = Series [ Sqrt [ 1 − x ] , { x , 0 , 20 } ] // Normal

1 − x 2 − x 2 8 − x 3 16 − 5 x 4 128 − 7 x 5 256 − 21 x 6 1024 − 33 x 7 2048 − 429 x 8 32768 − 715 x 9 65536 − 2431 x 10 262144 − 4199 x 11 524288 − 29393 x 12 4194304 − 52003 x 13 8388608 − 185725 x 14 33554432 − 334305 x 15 67108864 − 9694845 x 16 2147483648 − 17678835 x 17 4294967296 − 64822395 x 18 17179869184 − 119409675 x 19 34359738368 − 883631595 x 20 274877906944

In Fig. 3.61, observe that the series quickly converges to 1−x if x is small. He used this observation to use the series to approximate 3 by observing

3 = 4 ⋅ 3 4 = 2 1 − 1 4

and substituting x=1/4 into the first few terms of the series.

Figure 3.61 The Maclaurin series for 1−x converges quickly to 1−x when x is small (Virginia Polytechnic Institute & State University colors).

Plot [ { f [ x ] , Sqrt [ 1 − x ] } , { x , − 5 , 5 } , PlotRange → { − 5 , 5 } , AspectRatio → 1 ,

PlotStyle → { { Thickness [ . 01 ] , CMYKColor [ . 4 , 1 , . 5 , . 15 , . 5 ] } ,

{Thickness[.01],CMYKColor[0,.65,.9,0],Dashing[.02]}},AxesLabel→{x,y}]

Next, Newton looked at the circle x2−x+y2=(x−1/2)2+y2=(1/2)2 with center at (1/2,0) and radius 1/2. Refer to Fig. 3.62.

Figure 3.62 Newton's circle. The area of the green region is ∫01/4x−x2dx.

p1 = PolarPlot [ Cos [ r ] , { r , 0 , 2 Pi } , PlotStyle → { Thickness [ . 01 ] ,

CMYKColor [ . 4 , 1 , . 5 , . 15 ] } ] ;

p2 = Graphics [ { Thickness [ . 01 ] , CMYKColor [ 0 , . 65 , . 9 , 0 ] ,

Line [ { { 1 / 2 , 0 } , { 1 / 4 , Sqrt [ 3 ] / 4 } } ] , Line [ { { 1 / 4 , 0 } , { 1 / 4 , Sqrt [ 3 ] / 4 } } ] ,

CMYKColor [ 0 , . 18 , 1 , . 27 ] ,

Line [ { { 1 / 4 , − Sqrt [ 3 ] / 4 } , { 3 / 4 , Sqrt [ 3 ] / 4 } } ] ,

Line [ { { 1 / 2 , 0 } , { 3 / 4 , − Sqrt [ 3 ] / 4 } } ] } ] ;

p3 = Plot [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 1 / 4 } ,

PlotStyle → { Thickness [ . 01 ] , CMYKColor [ . 4 , 1 , . 5 , . 15 ] } ,

Filling → Bottom , FillingStyle → CMYKColor [ . 43 , 0 , . 14 , . 21 ] ] ;

Show [ p1 , p2 , p3 , AxesLabel → { x , y } , PlotLabel → “Newton's Circle” ]

Newton understood that the area of a circle with radius r is πr2 so it followed that the area of Newton's circle is π/4. Moreover, the area of the sector formed by the green region and the orange triangle is 1/6th the area of the circle so the area of the sector is Asector=π/24. The orange triangle has width 1/4 and height 3/4 so the orange triangle has area

A t r i a n g l e = 1 2 ⋅ 1 4 ⋅ 3 4 = 3 32 .

Because the top half of the circle is given by y=x−x2, the area of the green region, R, is AR=∫01/4x−x2dx. Newton then added and solved for π.

A s e c t o r = A R + A t r i a n g l e π 24 = ∫ 0 1 / 4 x − x 2 d x + 1 32 3 π = 24 ( ∫ 0 1 / 4 x − x 2 d x + 1 32 3 ) .

Newton already had an excellent approximation of 3 as described above. Thus, to obtain an excellent approximation of π he needed only to find an excellent approximation of ∫01/4x−x2dx. To obtain it, he started with his power series for 1−x and then observed that for 0⩽x⩽1, x−x2=xx−1=x1/2(1−x)1/2. With that in mind, he decided to take his power series for 1−x and multiply it by x1/2 to obtain a series that he predicted would converge to x−x2. Newton was right and it worked!

p1 = Series [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 9 } ] // Normal

x − x 3 / 2 2 − x 5 / 2 8 − x 7 / 2 16 − 5 x 9 / 2 128 − 7 x 11 / 2 256 − 21 x 13 / 2 1024 − 33 x 15 / 2 2048 − 429 x 17 / 2 32768

He then integrated his series (which was no longer a power series) term-by-term

Integrate [ p1 , x ] // Simplify

2 x 3 / 2 3 − x 5 / 2 5 − x 7 / 2 28 − x 9 / 2 72 − 5 x 11 / 2 704 − 7 x 13 / 2 1664 − 7 x 15 / 2 2560 − 33 x 17 / 2 17408 − 429 x 19 / 2 311296

and evaluated it if x=1/4 to obtain his approximation of ∫01/4x−x2dx.

There is no need to show the evaluation of the series if x=0 because the value of the series if x=0 is 0.

Integrate [ Series [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 9 } ] // Normal , x ] /. x → 1 / 4

1919013728784979 24995910798802944

We use N to obtain a numerical result.

Integrate [ Series [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 9 } ] // Normal , x ] /. x → 1 / 4 // N

0.0767731

We use Table to illustrate the steps Newton carried out by hand. In the first column, k, in the second column, the kth term of the series for x−x2 centered at x=0, in the third column, an antiderivative of the kth term of the series for x−x2 centered at x=0, and in the fourth column, the value of the antiderivative of the kth term of the series for x−x2 centered at x=0 evaluated at x=1/4.

Table [ { k , f [ k ] , x ∧ ( 1 / 2 ) f [ k ] , Integrate [ x ∧ ( 1 / 2 ) f [ k ] , x ] ,

Integrate [ x ∧ ( 1 / 2 ) f [ k ] , x ] /. x → 1 / 4 } , { k , 0 , 9 } ] // TableForm

0 1 x 2 x 3 / 2 3 1 12 1 − x 2 − x 3 / 2 2 − x 5 / 2 5 − 1 160 2 − x 2 8 − x 5 / 2 8 − x 7 / 2 28 − 1 3584 3 − x 3 16 − x 7 / 2 16 − x 9 / 2 72 − 1 36864 4 − 5 x 4 128 − 5 x 9 / 2 128 − 5 x 11 / 2 704 − 5 1441792 5 − 7 x 5 256 − 7 x 11 / 2 256 − 7 x 13 / 2 1664 − 7 13631488 6 − 21 x 6 1024 − 21 x 13 / 2 1024 − 7 x 15 / 2 2560 − 7 83886080 7 − 33 x 7 2048 − 33 x 15 / 2 2048 − 33 x 17 / 2 17408 − 33 2281701376 8 − 429 x 8 32768 − 429 x 17 / 2 32768 − 429 x 19 / 2 311296 − 429 163208757248 9 − 715 x 9 65536 − 715 x 19 / 2 65536 − 715 x 21 / 2 688128 − 715 1443109011456

Thus, Newton's approximation of π, π=24(∫01/4x−x2dx+1323) using his approximations of ∫01/4x−x2dx was given by

24 ( Sqrt [ 3 ] / 32 + Integrate [ Series [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 9 } ] // Normal , x ] /. x → 1 / 4 )

24 ( 1919013728784979 24995910798802944 + 3 32 )

We use N to approximate 3.

N [ 24 ( Sqrt [ 3 ] / 32 + Integrate [ Series [ Sqrt [ x − x ∧ 2 ] , { x , 0 , 9 } ] // Normal , x ] /. x → 1 / 4 ) ,

20 ]

3.1415926683631729578

When we compare Newton's approximation of π to our present day approximation of π, it is truly remarkable that Newton's approximation was accurate to seven decimal places (to the right of the decimal point).

N [ Pi , 20 ]

3.1415926535897932385

There are volumes of literature about the history of π, its applications, and how to efficiently and accurately approximate π (as well as other irrational and transcendental numbers).

Example 3.58

Determine whether the following series converge or diverge: (a) ∑k=0∞116k(48k+1−28k+4−18k+5−18k+6); (b) ∑k=0∞2k(k!)2(2k+1)!; (c) 229801∑k=0∞(4k)!(1103+26390k)(k!)43964k; (d) 12∑k=0∞(−1)k(6k)!(135910409+545140134k)(3k)!(k!)36403203k+3/2.

Solution

To determine whether each series converges or diverges, we use the Ratio Test. For (a), we see that limk→∞⁡|ak+1ak|=116. Because 1/16<1, the series converges (absolutely).

a [ k _ ] = 1 / 16 ∧ k ( 4 / ( 8 k + 1 ) − 2 / ( 8 k + 4 ) − 1 / ( 8 k + 5 ) − 1 / ( 8 k + 6 ) )

16 − k ( 4 1 + 8 k − 2 4 + 8 k − 1 5 + 8 k − 1 6 + 8 k )

s1 = FullSimplify [ a [ k + 1 ] / a [ k ] ]

( 1 + 2 k ) ( 3 + 4 k ) ( 1 + 8 k ) ( 5 + 8 k ) ( 318 + k ( 391 + 120 k ) ) 16 ( 3 + 2 k ) ( 7 + 4 k ) ( 9 + 8 k ) ( 13 + 8 k ) ( 47 + k ( 151 + 120 k ) )

Limit [ s1 , k → Infinity ]

1 16

In fact, in 1997, David Bailey, Peter Borwein, and Simon Plouffe (BBP) proved that

π = ∑ k = 0 ∞ 1 16 k ( 4 8 k + 1 − 2 8 k + 4 − 1 8 k + 5 − 1 8 k + 6 ) .

Notice that the second partial sum is accurate to three decimal places,

N[Sum[a[k],{k,0,1}],10]

3.141422466

and the seventh partial sum is accurate to ten decimal places.

N [ Sum [ a [ k ] , { k , 0 , 6 } ] , 15 ]

3.14159265357288

N [ Pi , 15 ]

3.14159265358979

Interestingly, in 2013, Fabrice Bellard improved on the BBP equation and proved that

π = 1 2 6 ∑ k = 0 ∞ ( − 1 ) k 2 10 k ( − 2 5 4 k + 1 − 1 4 k + 3 + 2 8 10 k + 1 − 2 6 10 k + 3 − 2 2 10 k + 5 − 2 2 10 k + 7 + 1 10 k + 9 ) .

(b) Using the Ratio Test, we see that the series converges because the limit of the ratio of successive terms is 1/2<1: limk→∞⁡|ak+1ak|=12.

a [ k _ ] = 2 ∧ k ( k ! ) ∧ 2 / ( 2 k + 1 ) !

2 k ( k ! ) 2 ( 1 + 2 k ) !

s1 = FullSimplify [ a [ k + 1 ] / a [ k ] ]

1 + k 3 + 2 k

Limit [ s1 , k → Infinity ]

1 2

Newton proved that this series converges to π/2:

π 2 = ∑ k = 0 ∞ 2 k ( k ! ) 2 ( 2 k + 1 ) ! .

Observe that the “patterns” of the terms of this series are much easier to see than in the example illustrating Newton's approximation of π using integral calculus, infinite series, and Euclidean geometry. You can use the 31st partial sum of this series to approximate π/2 to 10 decimal places.

N [ Sum [ a [ k ] , { k , 0 , 30 } ] , 10 ]

1.570796327

N [ Pi / 2 , 10 ]

1.570796327

(c) Using the Ratio Test, we see that

a [ k _ ] = 2 Sqrt [ 2 ] / 9801 ( 4 k ) ! ( 1103 + 26390 k ) / ( ( k ! ) ∧ 4396 ∧ ( 4 k ) )

2 3 2 − 8 k 99 − 2 − 4 k ( 1103 + 26390 k ) ( 4 k ) ! ( k ! ) 4

p1 = FullSimplify [ a [ k + 1 ] / a [ k ] ]

( 1 + 2 k ) ( 1 + 4 k ) ( 3 + 4 k ) ( 27493 + 26390 k ) 3073907232 ( 1 + k ) 3 ( 1103 + 26390 k )

Limit [ p1 , k → Infinity ]

1 96059601

limk→∞⁡|ak+1ak|=196059601<1 so the series converges absolutely. This series that was introduced by Srinivasa Ramanujan around 1910 converges to 1/π:

1 π = 2 2 9801 ∑ k = 0 ∞ ( 4 k ) ! ( 1103 + 26390 k ) ( k ! ) 4 396 4 k .

As with the other series in this example, the convergence is quick.

N [ Sum [ a [ k ] , { k , 0 , 2 } ] , 24 ]

0.318309886183790671537767

N [ 1 / Pi , 24 ]

0.318309886183790671537768

With just the third partial sum of the series, we see that it approximates 1/π to 23 decimal places.

(d) Using the Ratio Test, we see that the series converges. This formulation of π was introduced by brothers David Chudnosky and Gregory Chudnosky in 1989:

1 π = 12 ∑ k = 0 ∞ ( − 1 ) k ( 6 k ) ! ( 135910409 + 545140134 k ) ( 3 k ) ! ( k ! ) 3 640320 3 k + 3 / 2 .

a [ k _ ] = 12 ( − 1 ) ∧ k ( 6 k ) ! ( 13591409 + 545140134 k ) / ( ( 3 k ) ! ( k ! ) ∧ 3640320 ∧ ( 3 k + 3 / 2 ) )

( − 1 ) k 3 − 1 2 − 3 k 4 − 7 2 − 9 k 3335 − 3 2 − 3 k ( 13591409 + 545140134 k ) ( 6 k ) ! ( k ! ) 3 ( 3 k ) !

p1 = Abs [ FullSimplify [ a [ k + 1 ] / a [ k ] ] ]

Abs [ ( 1 + 2 k ) ( 1 + 6 k ) ( 5 + 6 k ) ( 558731543 + 545140134 k ) ( 1 + k ) 3 ( 13591409 + 545140134 k ) ] 10939058860032000

Limit[p1,k→Infinity]

1 151931373056000

The series converges very quickly. Using the third partial sum, we see that it yields an approximation of 1/π accurate to 42 decimal places.

N [ Sum [ a [ k ] , { k , 0 , 2 } ] , 42 ]

0.318309886183790671537767526745028724068919

N [ 1 / Pi , 42 ]

0.318309886183790671537767526745028724068919

To use the Chudnosky's series to approximate π, we first rewrite it as

1 π = 1 42680 10005 ∑ k = 0 ∞ ( − 1 ) k ( 6 k ) ! ( 13591409 + 545140134 k ) ( 3 k ) ! ( k ! ) 3 640320 3 k .

(3.33)

Now let A=∑k=0∞ak=∑k=0∞(−1)k(6k)!(3k)!(k!)36403203k and B=∑k=0∞bk=∑k=0∞(−1)k(6k)!k(3k)!(k!)36403203k. Then (3.33) becomes

1 π = 1351409 A + 545140134 B 42680 10005

(3.34)

and reciprocating gives us

π = 42680 10005 1351409 A + 545140134 B .

(3.35)

Observe that ak=(−1)k(6k)!(3k)!(k!)36403203k, bk=k⋅ak, and

a k a k − 1 = − ( 6 k − 1 ) ( 6 k − 5 ) ( 2 k − 1 ) 10939058860032000 k 2 a k = − ( 6 k − 1 ) ( 6 k − 5 ) ( 2 k − 1 ) 10939058860032000 k 2 a k − 1 .

Then, if we let An=∑k=0nak,

A n = ∑ k = 0 n − 1 a k − ( 6 n − 1 ) ( 6 n − 5 ) ( 2 n − 1 ) 10939058860032000 ( n − 1 ) n 2 a n − 1 = A n − 1 − ( 6 n − 1 ) ( 6 n − 5 ) ( 2 n − 1 ) 10939058860032000 ( n − 1 ) n 2 a n − 1

and if we let Bn=∑k=0nbk, Bn=∑k=0nkak=∑k=0n−1kak+nan=Bn−1+nan. Then, we approximate π with π≈=42680100051351409An+545140134Bn.

Remove [ a , b , pi , capa , capb ]

a [ 0 ] = 1 ;

a [ n _ ] := a [ n ] = N [ − ( − 1 + 2 n ) ( − 5 + 6 n ) ( − 1 + 6 n ) 10939058860032000 n 3 a [ n − 1 ] , 250 ] ;

b[n_]:=b[n]=na[n]

capa [ 0 ] = a [ 0 ] ;

capb [ 0 ] = b [ 0 ] ;

capa [ n _ ] := capa [ n ] = capa [ n − 1 ] + a [ n ]

capb [ n _ ] := capb [ n ] = capb [ n − 1 ] + b [ n ]

pi [ n _ ] := N [ 426880 Sqrt [ 10005 ] / ( 13591409 capa [ n ] + 545140134 capb [ n ] ) , 250 ]

We illustrate the quick convergence using A10 and B10. Observe that the approximation gives us an approximation of π accurate to at least 155 digits.

N [ pi [ 10 ] , 155 ]

3.1415926535897932384626433832795028841971693993751058209749 445923078164062862089986280348253421170679821480865132823066 470938446095505822317253594081284811

N [ Pi , 155 ]

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811 □

Closely related to these series is the Madhava–Leibniz series, π=12∑k=0∞(−1)k2k+13k. It is a good exercise to verify that the first 21 terms of the Madhava–Leibniz series compute an approximation of π accurate to eleven decimal places.

3.4.8 Other Series

In calculus, we learn that the power series f(x)=∑k=0∞ak(x−x0)k is differentiable and integrable on its interval of convergence. However, for series that are not power series this result is not generally true. For example, in more advanced courses, we learn that the function

f ( x ) = ∑ k = 0 ∞ 1 2 k sin ⁡ ( 3 k x )

is continuous for all values of x but nowhere differentiable. We can use Mathematica to help us see why this function is not differentiable. Let

f n ( x ) = ∑ k = 0 n 1 2 k sin ⁡ ( 3 k x ) .

Notice that fn(x) is defined recursively by f0(x)=sin⁡x and fn(x)=fn−1(x)+12nsin⁡(3nx). We use Mathematica to recursively define fn(x).

Clear [ f ]

f [ 0 ] = Sin [ x ] ;

f[k_]:=f[k]=f[k−1]+Sin[3kx]2k

We define fn(x) using the form

f[n_]:=f[n]=...

so that Mathematica “remembers” the values it computes. Thus, to compute f[5], Mathematica uses the previously computed values, namely f[4], to compute f[5]. Note that we can produce the same results by defining fn(x) with the command

f[n_]:=...

However, the disadvantage of defining fn(x) in this manner is that Mathematica does not “remember” the previously computed values and thus takes longer (and uses more memory) to compute fn(x) for larger values of n.

Next, we use Table to generate f3(x), f6(x), f9(x), and f12(x).

tograph = Table [ f [ n ] , { n , 3 , 12 , 3 } ]

We now graph each of these functions and show the results as a graphics array with GraphicsGrid in Fig. 3.63. (Note that you do not need to include the option DisplayFunction->Identity to suppress the resulting output unless you forget to include the semi-colon at the end of the command.)

Figure 3.63 Approximating a function that is continuous everywhere but nowhere differentiable (Duke University colors).

graphs = Table [ Plot [ Evaluate [ tograph [ [ i ] ] ] , { x , 0 , 3 π } ,

PlotStyle → CMYKColor [ 1 , . 75 , . 06 , . 24 ] ] , { i , 1 , 4 } ] ;

toshow = Partition [ graphs , 2 ] ;

Show [ GraphicsGrid [ toshow ] ]

From these graphs, we see that for large values of n, the graph of fn(x), although actually smooth, appears “jagged” and thus we might suspect that f(x)=limn→∞ fn(x)=∑k=0∞12ksin⁡(3kx) is indeed continuous everywhere but nowhere differentiable.

3.5 Multi-Variable Calculus

Mathematica is useful in investigating functions involving more than one variable. In particular, the graphical analysis of functions that depend on two (or more) variables is enhanced with the help of Mathematica's graphics capabilities.

3.5.1 Limits of Functions of Two Variables

Mathematica's graphics and numerical capabilities are helpful in investigating limits of functions of two variables.

Example 3.59

Show that the limit lim(x,y)→(0,0)⁡x2−y2x2+y2 does not exist.

Solution

We begin by defining f(x,y)=(x2−y2)/(x2+y2). Next, we use Plot3D to graph z=f(x,y) for −1/2⩽x⩽1/2 and −1/2⩽y⩽1/2. ContourPlot is used to graph several level curves on the same rectangle. (See Fig. 3.64.) (To define a function of two variables, f(x,y)=expressioninxandy, enter f[x_,y_]=expression in x and y. Plot3D[f[x,y],{a,x,b},{y,c,d}] generates a basic graph of z=f(x,y) for a⩽x⩽b and c⩽y⩽d)

Figure 3.64 (a) Three-dimensional and (b) contour plots of f(x,y) (North Carolina State University colors).

Clear [ f ]

f [ x _ , y _ ] = ( x ∧ 2 − y ∧ 2 ) / ( x ∧ 2 + y ∧ 2 ) ;

p1 = Plot3D [ f [ x , y ] , { x , − 2 , 2 } , { y , − 2 , 2 } , PlotPoints → 40 ,

PlotStyle → { Opacity [ . 8 ] , CMYKColor [ 0 , 1 , . 81 , . 04 ] } ,

MeshFunctions → { # 3 & } , MeshStyle → { Black , Thickness [ . 005 ] } ] ;

c1 = ContourPlot [ f [ x , y ] , { x , − 2 , 2 } , { y , − 2 , 2 } , ContourShading → False ,

Axes → Automatic , AxesOrigin → { 0 , 0 } , PlotPoints → 60 , Contours → 20 ,

ContourStyle → CMYKColor [ 0 , 1 , . 81 , . 04 ] , Frame → False ,

AxesLabel → { x , y } ] ;

Show [ GraphicsRow [ { p1 , c1 } ] ]

From the graph of the level curves, we suspect that the limit does not exist because we see that near (0,0), z=f(x,y) attains many different values. We obtain further evidence that the limit does not exist by computing the value of z=f(x,y) for various points chosen randomly near (0,0). We use Table and RandomReal to generate 10 ordered pairs (x,y) for x and y “close to” 0. Because RandomReal is included in the calculation, your results will almost certainly be different from those here.

When you slide the cursor over the contours in the contour plot, the contour values are displayed.

pts = Table [ RandomReal [ { − 10 − i , 10 − i } ] , { i , 1 , 10 } , { 2 } ]

{ { − 0.0240696 , 0.0964425 } , { 0.00929171 , 0.00582256 } , { − 0.000875533 , − 0.000861975 } , { 4.93995041959814 ` * − ∧ 6 , 0.0000585256 } , { − 5.7785374509223095 ` * − ∧ 6 , 6.7869167280744185 ` * − ∧ 6 } , { − 6.964772344758519 ` * − ∧ 7 , − 9.809015465572932 ` * − ∧ 8 } , { 1.669922073323435 ` * − ∧ 8 , − 7.416883164217774 ` * − ∧ 8 } , { 4.256718340841544 ` * − ∧ 9 , 2.2211732111688833 ` * − ∧ 10 } , { 4.8322758801790404 ` * − ∧ 11 , − 4.063639702523537 ` * − ∧ 10 } , { 2.96738846034934 ` * − ∧ 11 , − 9.412273092884795 ` * − ∧ 11 } }

Next, we define a function g that given an ordered pair (x,y) ({x,y} in Mathematica), g((x,y)) returns the ordered triple (x,y,f(x,y)) ({x,y,f[x,y]} in Mathematica).

g [ { x _ , y _ } ] = { x , y , f [ x , y ] }

{ x , y , x 2 − y 2 x 2 + y 2 }

We then use Map to apply g to the list pts.

Map [ g , pts ] // TableForm

−0.02406960.0964425−0.8827290.009291710.005822560.436083−0.000875533−0.0008619750.01560524.93995041959814`*−∧60.0000585256−0.985852−5.7785374509223095`*−∧66.7869167280744185`*−∧6−0.159473−6.964772344758519`*−∧7−9.809015465572932`*−∧80.9611011.669922073323435`*−∧8−7.416883164217774`*−∧8−0.9035054.256718340841544`*−∧92.2211732111688833`*−∧100.9945694.8322758801790404`*−∧11−4.063639702523537`*−∧10−0.9721132.96738846034934`*−∧11−9.412273092884795`*−∧11−0.819184

From the third column, we see that z=f(x,y) does not appear to approach any particular value for points chosen randomly near (0,0). In fact, along the line y=mx we see that f(x,y)=f(x,mx)=(1−m2)/(1+m2). Hence as (x,y)→(0,0) along y=mx, f(x,y)=f(x,mx)→1−m21+m2. Thus, f(x,y) does not have a limit as (x,y)→(0,0).

We choose lines of the form y=mx because near (0,0) the level curves of z=f(x,y) look like lines of the form y=mx.

v1 = Simplify [ f [ x , m x ] ]

1 − m 2 1 + m 2

v1 /. m → 0

v1 /. m → 1

v1 /. m → 1 / 2

1

0

35 □

In some cases, you can establish that a limit does not exist by converting to polar coordinates. For example, in polar coordinates, f(x,y)=x2−y2x2+y2 becomes f(rcos⁡θ,rsin⁡θ)=2cos2⁡θ−1

Simplify [ f [ r Cos [ t ] , r Sin [ t ] ] ]

Cos [ 2 t ]

and

lim ( x , y ) → ( 0 , 0 ) ⁡ f ( x , y ) = lim r → 0 ⁡ f ( r cos ⁡ θ , r sin ⁡ θ ) = lim r → 0 ⁡ 2 cos 2 ⁡ θ − 1 = 2 cos 2 ⁡ θ − 1 = cos ⁡ 2 θ

depends on θ.

3.5.2 Partial and Directional Derivatives

Partial derivatives of functions of two or more variables are computed with Mathematica using D. For z=f(x,y),

1.  D[f[x,y],x] computes ∂f∂x=fx(x,y),

2.  D[f[x,y],y] computes ∂f∂y=fy(x,y),

3.  D[f[x,y],{x,n}] computes ∂nf∂xn,

4.  D[f[x,y],y,x] computes ∂2f∂y∂x=fxy(x,y), and

5.  D[f[x,y],{x,n},{y,m}] computes ∂n+mf∂nx∂my.

6.  You can use the button located on the BasicMathInput palette to create templates to compute partial derivatives.

The calculations are carried out similarly for functions of more than two variables.

Example 3.60

Calculate fx(x,y), fy(x,y), fxy(x,y), fyx(x,y), fxx(x,y), and fyy(x,y) if f(x,y)=sin⁡x2+y2+1.

Solution

After defining f(x,y)=sin⁡x2+y2+1,

f [ x _ , y _ ] = Sin [ Sqrt [ x ∧ 2 + y ∧ 2 + 1 ] ] ;

we illustrate the use of D to compute the partial derivatives. Entering

D [ f [ x , y ] , x ]

x Cos [ 1 + x 2 + y 2 ] 1 + x 2 + y 2

computes fx(x,y). Entering

D [ f [ x , y ] , y ]

y Cos [ 1 + x 2 + y 2 ] 1 + x 2 + y 2

computes fy(x,y). Entering

D [ f [ x , y ] , x , y ] // Together

− x y ( Cos [ 1 + x 2 + y 2 ] + 1 + x 2 + y 2 Sin [ 1 + x 2 + y 2 ] ) ( 1 + x 2 + y 2 ) 3 / 2

computes fyx(x,y). Entering

D [ f [ x , y ] , y , x ] // Together

− x y ( Cos [ 1 + x 2 + y 2 ] + 1 + x 2 + y 2 Sin [ 1 + x 2 + y 2 ] ) ( 1 + x 2 + y 2 ) 3 / 2

computes fxy(x,y). Remember that under appropriate assumptions, fxy(x,y)=fyx(x,y). Entering

D [ f [ x , y ] , { x , 2 } ] // Together

Cos [ 1 + x 2 + y 2 ] + y 2 Cos [ 1 + x 2 + y 2 ] − x 2 1 + x 2 + y 2 Sin [ 1 + x 2 + y 2 ] ( 1 + x 2 + y 2 ) 3 / 2

computes fxx(x,y). Entering

D [ f [ x , y ] , { y , 2 } ] // Together

Cos [ 1 + x 2 + y 2 ] + x 2 Cos [ 1 + x 2 + y 2 ] − y 2 1 + x 2 + y 2 Sin [ 1 + x 2 + y 2 ] ( 1 + x 2 + y 2 ) 3 / 2

computes fyy(x,y).  □

The directional derivative of z=f(x,y) in the direction of the unit vector u=cos⁡θi+sin⁡θj is

D u f ( x , y ) = f x ( x , y ) cos ⁡ θ + f y ( x , y ) sin ⁡ θ ,

provided that fx(x,y) and fy(x,y) both exist.

The vectors i and j are defined by i=〈1,0〉 and j=〈0,1〉.

If fx(x,y) and fy(x,y) both exist, the gradient of f(x,y) is the vector-valued function

▽ f ( x , y ) = f x ( x , y ) i + f y ( x , y ) j = 〈 f x ( x , y ) , f y ( x , y ) 〉 .

Notice that if u=〈cos⁡θ,sin⁡θ〉,

D u f ( x , y ) = ▽ f ( x , y ) ⋅ 〈 cos ⁡ θ , sin ⁡ θ 〉 .

Let (x0,y0) be a point in the domain of the differentiable function z=f(x,y). At (x0,y0), f(x,y) increases most rapidly in the direction of the gradient, ▽f(x,y), evaluated at (x0,y0). Similarly, f(x,y) decreases most rapidly in the direction of −▽f(x,y). Consequently, the gradient is perpendicular to the level curves of f(x,y).

Calculus of vector-valued functions is discussed in more detail in Chapter 5.

StreamPlot[{f(x,y),g(x,y)},{x,a,b},{y,c,d}]

generates a stream plot of the vector field <f(x,y),g(x,y)>. Thus,

StreamPlot[{D[f[x,y],x],D[f[x,y],y]},{x,a,b},{y,c,d}]

generates a stream plot of the gradient of f(x,y).

Vectors and vector-valued functions will be discussed in more detail in Chapters 4 and 5.

Example 3.61

Let f(x,y)=6x2y−3x4−2y3. (a) Find Duf(x,y) in the direction of v=〈3,4〉. (b) Compute D〈3/5,4/5〉f(139+33,1). (c) Find an equation of the line tangent to the graph of 6x2y−3x4−2y3=0 at the point (139+33,1).

Solution

After defining f(x,y)=6x2y−3x4−2y3, we graph z=f(x,y) with Plot3D in Fig. 3.65, illustrating the PlotPoints, PlotRange, and ViewPoint options.

Figure 3.65 f(x,y)=6x²y − 3x⁴ − 2y³ for −2 ⩽ x ⩽ 2 and −2 ⩽ y ⩽ 3 (Rice University colors).

f [ x _ , y _ ] = 6 x ∧ 2 y − 3 x ∧ 4 − 2 y ∧ 3 ;

Plot3D [ f [ x , y ] , { x , − 2 , 2 } , { y , − 2 , 3 } , PlotPoints → 50 ,

PlotRange → { { − 2 , 2 } , { − 2 , 3 } , { − 2 , 2 } } ,

PlotStyle → { CMYKColor [ 1 , . 93 , . 28 , . 22 ] , Opacity [ . 5 ] } ,

BoxRatios → { 1 , 1 , 1 } , ViewPoint → { 1.887 , 2.309 , 1.6 } ,

ClippingStyle → None ]

(a) A unit vector, u, in the same direction as v is

u = 〈 3 3 2 + 4 2 , 4 3 2 + 4 2 〉 = 〈 3 5 , 4 5 〉 .

v = { 3 , 4 } ;

u = v / Sqrt [ v . v ]

{ 3 5 , 4 5 }

Then, Duf(x,y)=〈fx(x,y),fy(x,y)〉⋅u, calculated in du.

gradf = { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] }

{ − 12 x 3 + 12 x y , 6 x 2 − 6 y 2 }

du = Simplify [ gradf . u ]

− 12 5 ( − 2 x 2 + 3 x 3 − 3 x y + 2 y 2 )

(b) D〈3/5,4/5〉f(139+33,1) is calculated by evaluating du if x=139+33 and y=1.

du1 = du /. { x → 1 / 3 Sqrt [ 9 + 3 Sqrt [ 3 ] ] , y → 1 } // Simplify

4 5 ( 2 3 − 3 3 + 3 )

(c) The gradient is evaluated if x=139+33 and y=1.

nvec = gradf /. { x → 1 / 3 Sqrt [ 9 + 3 Sqrt [ 3 ] ] , y → 1 } // Simplify

{−43+3,23}

Generally, ▽f(x,y) is perpendicular to the level curves of z=f(x,y), so

nvec = ▽ f ( 1 3 9 + 3 3 , 1 ) = 〈 f x ( 1 3 9 + 3 3 , 1 ) , f y ( 1 3 9 + 3 3 , 1 ) 〉

is perpendicular to f(x,y)=0 at the point (139+33,1). Thus, an equation of the line tangent to the graph of f(x,y)=0 at the point (139+33,1) is

f x ( 1 3 9 + 3 3 , 1 ) ( x − 1 3 9 + 3 3 ) + f y ( 1 3 9 + 3 3 , 1 ) ( y − 1 ) = 0 ,

which we solve for y with Solve. We confirm this result by graphing f(x,y)=0 using ContourPlot in conf and then graphing the tangent line in tanplot. tanplot and conf are shown together with Show in Fig. 3.66.

An equation of the line L containing (x0,y0) and perpendicular to n=〈a,b〉 is a(x−x0)+b(y−y0)=0.

Figure 3.66 Level curve of f(x,y) together with a tangent line.

conf = ContourPlot [ f [ x , y ] = = 0 , { x , − 2 , 2 } , { y , − 2 , 2 } , PlotPoints → 60 ,

ContourShading → False , Frame → False , Axes → Automatic ,

ContourStyle - > CMYKColor [ 1 , . 93 , . 28 , . 22 ] ,

AxesOrigin → { 0 , 0 } ] ;

tanline=Solve[nvec[[1]](x−1/3Sqrt[9+3Sqrt[3]])+nvec[[2]](y−1)==0,y]

{ { y → − 2 + 3 − 2 3 + 3 x 3 } }

tanplot = Plot [ Evaluate [ y /. tanline ] , { x , − 2 , 2 } ,

PlotStyle - > CMYKColor [ 0 , 0 , 0 , . 77 ] ] ;

Show [ conf , tanplot , PlotRange → { { − 2 , 2 } , { − 2 , 3 } } , AspectRatio → Automatic ]

More generally, we use ContourPlot together with the StreamPlot function to illustrate that the gradient vectors are perpendicular to the level curves of z=f(x,y) in Fig. 3.67.

Figure 3.67 (a) Level curves of z = f(x,y). (b) Gradient field of z = f(x,y). (c) The gradient together with several level curves.

p1 = ContourPlot [ f [ x , y ] , { x , − 2 , 2 } , { y , − 2 , 3 } ,

Contours → 25 , ContourStyle → CMYKColor [ 0 , 0 , 0 , . 77 ] ,

ContourShading → False , Frame → False , Axes → Automatic ,

AxesOrigin → { 0 , 0 } , AxesLabel → { x , y } ,

PlotLabel → “(a)” ]

p2 = StreamPlot [ Evaluate [ { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] } ] ,

{ x , − 2 , 2 } , { y , − 2 , 3 } , StreamStyle → CMYKColor [ 1 , . 93 , . 28 , . 22 ] ,

Frame → False , Axes → Automatic , AxesLabel → { x , y } ,

PlotLabel → “(b)” ]

p3 = Show [ p1 , p2 , PlotLabel → “(c)” ]

Show [ GraphicsRow [ { p1 , p2 , p3 } ] ]

In Figs. 3.67 (b) and (c), observe that if the gradient vectors are all converging to a point, ▽f(x,y)=0 and the point will be a relative maximum. At (0,0), ▽f(x,y)=0 but at this point, some gradient vectors are pointing towards the origin and some away so (0,0) is a saddle. Compare to Fig. 3.65. Note that if ▽f(x,y)=0 and all the gradient vectors were pointing away from the point, the point would be a relative minimum. □

Example 3.62

Let

f ( x , y ) = ( y − 1 ) 2 e − ( x + 1 ) 2 − y 2 − 10 3 ( − x 5 + 1 5 y − y 3 ) e − x 2 − y 2 − 1 9 e − x 2 − ( y + 1 ) 2 .

Calculate ▽f(x,y) and then graph ▽f(x,y) together with several level curves of f(x,y).

Solution

We begin by defining and graphing z=f(x,y) with Plot3D in Fig. 3.68 (a).

Figure 3.68 (a) f(x,y) for −3 ⩽ x ⩽ 3 and −3 ⩽ y ⩽ 2. (b) Contour plot of f(x,y) along with several gradient vectors (University of Washington colors).

Clear [ f ]

f [ x _ , y _ ] = ( y − 1 ) ∧ 2 Exp [ − ( x + 1 ) ∧ 2 − y ∧ 2 ] −

10 / 3 ( − x ∧ 5 + 1 / 5 y − y ∧ 3 ) Exp [ − x ∧ 2 − y ∧ 2 ] −

1 / 9 Exp [ − x ∧ 2 − ( y + 1 ) ∧ 2 ] ;

p1 = Plot3D [ f [ x , y ] , { x , − 3 , 3 } , { y , − 3 , 3 } ,

PlotStyle → { CMYKColor [ 0 , . 13 , . 43 , . 13 ] , Opacity [ . 5 ] } ,

ViewPoint → { − 1.99 , 2.033 , 1.833 } ,

PlotRange → All , PlotLabel → “(a)” ] ;

conf = ContourPlot [ f [ x , y ] , { x , − 3 , 3 } , { y , − 3 , 3 } ,

PlotPoints → 60 , Contours → 30 , ContourShading → False ,

Frame → False , Axes → Automatic , ContourStyle → CMYKColor [ . 93 , 1 , . 18 , . 21 ] ,

AxesOrigin → { 0 , 0 } ] ;

In the three-dimensional plot, notice that z appears to have six relative extrema: three relative maxima and three relative minima. We also graph several level curves of f(x,y) with ContourPlot and name the resulting graphic conf.

Next we calculate fx(x,y) and fy(x,y) using Simplify and D. The gradient is the vector-valued function 〈fx(x,y),fy(x,y)〉.

gradf = { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] } // Simplify

{ − 2 9 e − 2 x − x 2 − ( 1 + y ) 2 ( − e 2 x x + 9 e 2 y ( 1 + x ) ( − 1 + y ) 2 + 3 e 1 + 2 x + 2 y x ( − 25 x 3 + 10 x 5 − 2 y + 10 y 3 ) ) , − 2 9 e − 2 x − x 2 − ( 1 + y ) 2 ( − e 2 x ( 1 + y ) + 9 e 2 y ( 1 − 2 y 2 + y 3 ) + e 1 + 2 x + 2 y ( 3 + 30 x 5 y − 51 y 2 + 30 y 4 ) ) }

We use StreamPlot to graph the gradient naming the resulting graphic gradfplot. gradfplot and conf are displayed together using Show in Fig. 3.68 (b).

gradfplot = StreamPlot [ gradf , { x , − 3 , 3 } , { y , − 3 , 3 } ,

StreamStyle → CMYKColor [ . 3 , . 35 , . 6 , 0 ]

];

Show [ GraphicsRow [ { p1 , Show [ conf , gradfplot , PlotLabel → “(b)” ] } ] ]

In the result (see Fig. 3.68 (b)), notice that the gradient is perpendicular to the level curves; the gradient is pointing in the direction of maximal increase of z=f(x,y). □

Classifying Critical Points

Let z=f(x,y) be a real-valued function of two variables with continuous second-order partial derivatives. A critical point of z=f(x,y) is a point (x0,y0) in the interior of the domain of z=f(x,y) for which

f x ( x 0 , y 0 ) = 0 and f y ( x 0 , y 0 ) = 0 ,

or, equivalently, ▽f(x0,y0)=0. Critical points are classified by the Second Derivatives (or Partials) test.

Theorem 3.18

Second Derivatives Test

Let (x0,y0) be a critical point of a function z=f(x,y) of two variables and let

d = f x x ( x 0 , y 0 ) f y y ( x 0 , y 0 ) − [ f x y ( x 0 , y 0 ) ] 2 .

(3.36)

1.  If d>0 and fxx(x0,y0)>0, then z=f(x,y) has a relative (or local) minimum at (x0,y0).

2.  If d>0 and fxx(x0,y0)<0, then z=f(x,y) has a relative (or local) maximum at (x0,y0).

3.  If d<0, then z=f(x,y) has a saddle point at (x0,y0).

4.  If d=0, no conclusion can be drawn and (x0,y0) is called a degenerate critical point.

Example 3.63

Find the relative maximum, relative minimum, and saddle points of f(x,y)=−2x2+x4+3y−y3.

Solution

After defining f(x,y), the critical points are found with Solve and named critpts.

f [ x _ , y _ ] = − 2 x ∧ 2 + x ∧ 4 + 3 y − y ∧ 3 ;

critpts = Solve [ { D [ f [ x , y ] , x ] = = 0 , D [ f [ x , y ] , y ] = = 0 } , { x , y } ]

{ { x → − 1 , y → − 1 } , { x → 0 , y → − 1 } , { x → 1 , y → − 1 } , { x → − 1 , y → 1 } , { x → 0 , y → 1 } , { x → 1 , y → 1 } }

We then define dfxx. Given (x0,y0), dfxx(x0,y0) returns the ordered quadruple x0, y0, equation (3.36) evaluated at (x0,y0), and fxx(x0,y0).

dfxx [ x0 _ , y0 _ ] = { x0 , y0 ,

D [ f [ x , y ] , { x , 2 } ] D [ f [ x , y ] , { y , 2 } ] − D [ f [ x , y ] , x , y ] ∧ 2 /.

{ x → x0 , y → y0 } , D [ f [ x , y ] , { x , 2 } ] /. { x → x0 , y → y0 } }

{ x 0 , y 0 , − 6 ( − 4 + 12 x 0 2 ) y 0 , − 4 + 12 x 0 2 }

For example,

dfxx [ 0 , 1 ]

{ 0 , 1 , 24 , − 4 }

shows us that a relative maximum occurs at (0,1). We then use /. (ReplaceAll) to substitute the values in each element of critpts into dfxx.

dfxx [ x , y ] /. critpts // TableForm

− 1 − 1 48 8 0 − 1 − 24 − 4 1 − 1 48 8 − 1 1 − 48 8 0 1 24 − 4 1 1 − 48 8

From the result, we see that (0,1) results in a relative maximum, (0,−1) results in a saddle, (1,1) results in a saddle, (1,−1) results in a relative minimum, (−1,1) results in a saddle, and (−1,−1) results in a relative minimum. We confirm these results graphically with a three-dimensional plot generated with Plot3D and a contour plot together with the gradient generated with ContourPlot and StreamPlot in Fig. 3.69.

Figure 3.69 (a) Three-dimensional and (b) contour plots of f(x,y) (University of Florida colors).

p1 = Plot3D [ f [ x , y ] , { x , − 3 / 2 , 3 / 2 } , { y , − 3 / 2 , 3 / 2 } , PlotPoints → 40 ,

PlotStyle → { CMYKColor [ 0 , . 7 , 1 , 0 ] , Opacity [ . 4 ] } ,

PlotLabel → “(a)” ] ;

p2 = ContourPlot [ f [ x , y ] , { x , − 3 / 2 , 3 / 2 } , { y , − 3 / 2 , 3 / 2 } ,

PlotPoints → 40 , ContourShading → False , Contours → 20 ,

Frame → False , Axes → Automatic , AxesOrigin → { 0 , 0 } ,

ContourStyle → CMYKColor [ 1 , . 6 , 0 , . 2 ] ,

AxesLabel → { x , y } ] ;

p3 = StreamPlot [ Evaluate [ { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] } ] , { x , − 3 / 2 , 3 / 2 } ,

{ y , − 3 / 2 , 3 / 2 } ,

StreamStyle → CMYKColor [ . 1 , . 2 , . 4 , . 1 ] ] ;

Show [ GraphicsRow [ { p1 , Show [ p2 , p3 , PlotLabel → “(b)” ] } ] ]

In the contour plot, notice that near relative extrema, the level curves look like circles while near saddles they look like hyperbolas. Further, at the maximums, the gradient is directed towards the point while at the minimums it is directed away from the point. □

If the Second Derivatives Test fails, graphical analysis is especially useful.

Example 3.64

Find the relative extrema and saddle points of f(x,y)=x2+x2y2+y4.

Solution

Initially we proceed in the exact same manner as in the previous example: we define f(x,y) and compute the critical points. Several complex solutions are returned, which we ignore.

f [ x _ , y _ ] = x ∧ 2 + x ∧ 2 y ∧ 2 + y ∧ 4 ;

critpts = Solve [ { D [ f [ x , y ] , x ] = = 0 , D [ f [ x , y ] , y ] = = 0 } , { x , y } ]

{{x→0,y→0},{x→0,y→0},{x→−2,y→−i},{x→2,y→−i},{x→−2,y→i},{x→2,y→i}}

We then compute the value of (3.36) at the real critical point, and the value of fxx(x,y) at this critical point.

dfxx [ x0 _ , y0 _ ] = { x0 , y0 ,

D [ f [ x , y ] , { x , 2 } ] D [ f [ x , y ] , { y , 2 } ] − D [ f [ x , y ] , x , y ] ∧ 2 /.

{ x → x0 , y → y0 } , D [ f [ x , y ] , { x , 2 } ] /. { x → x0 , y → y0 } }

{ x 0 , y 0 , − 16 x 0 2 y 0 2 + ( 2 + 2 y 0 2 ) ( 2 x 0 2 + 12 y 0 2 ) , 2 + 2 y 0 2 }

dfxx [ 0 , 0 ]

{ 0 , 0 , 0 , 2 }

The result shows us that the Second Derivatives Test fails at (0,0).

p1 = Plot3D [ f [ x , y ] , { x , − 1 , 1 } , { y , − 1 , 1 } , BoxRatios → Automatic ,

PlotStyle → { CMYKColor [ 1 , . 56 , 0 , . 34 ] , Opacity [ . 3 ] } ,

MeshFunctions - > { # 3 & } , MeshStyle → { CMYKColor [ 0 , . 18 , 1 , . 15 ] , Thickness [ . 0075 ] } ,

PlotLabel → “(a)” ] ;

p2 = ContourPlot [ f [ x , y ] , { x , − 1 , 1 } , { y , − 1 , 1 } , PlotPoints → 40 ,

Contours → 20 , ContourShading → False , Frame → False ,

Axes → Automatic , AxesOrigin → { 0 , 0 } , ContourStyle →

{ { CMYKColor [ 0 , . 18 , 1 , . 15 ] , Thickness [ . 0075 ] } } ] ;

p3 = StreamPlot [ Evaluate [ { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] } ] , { x , − 1 , 1 } ,

{ y , − 1 , 1 } ,

StreamStyle → CMYKColor [ . 03 , 1 , . 7 , . 12 ] ] ;

Show [ GraphicsRow [ { p1 , Show [ p2 , p3 , PlotLabel → b ] } ] ]

However, the contour plot and the stream plot of the gradient of f(x,y) near (0,0) indicates that an extreme value occurs at (0,0). Because the gradient vectors are all pointing away from the origin, (0,0) is a relative minimum which is confirmed in the three-dimensional plot. It is also an absolute minimum. (See Fig. 3.70.) □

Figure 3.70 (a) Three-dimensional and (b) contour plots of f(x,y) (University of California-Davis colors).

Tangent Planes

Let z=f(x,y) be a real-valued function of two variables. If both fx(x0,y0) and fy(x0,y0) exist, then an equation of the plane tangent to the graph of z=f(x,y) at the point (x0,y0,f(x0,y0)) is given by

f x ( x 0 , y 0 ) ( x − x 0 ) + f y ( x 0 , y 0 ) ( y − y 0 ) − ( z − z 0 ) = 0 ,

(3.37)

where z0=f(x0,y0). Solving for z yields the function (of two variables)

z = f x ( x 0 , y 0 ) ( x − x 0 ) + f y ( x 0 , y 0 ) ( y − y 0 ) + z 0 .

(3.38)

Symmetric equations of the line perpendicular to the surface z=f(x,y) at the point (x0,y0,z0) are given by

x − x 0 f x ( x 0 , y 0 ) = y − y 0 f y ( x 0 , y 0 ) = z − z 0 − 1

(3.39)

and parametric equations are

{ x = x 0 + f x ( x 0 , y 0 ) t y = y 0 + f y ( x 0 , y 0 ) t z = z 0 − t .

(3.40)

The plane tangent to the graph of z=f(x,y) at the point (x0,y0,f(x0,y0)) is the “best” linear approximation of z=f(x,y) near (x,y)=(x0,y0) in the same way as the line tangent to the graph of y=f(x) at the point (x0,f(x0)) is the “best” linear approximation of y=f(x) near x=x0.

Example 3.65

Find an equation of the plane tangent and normal line to the graph of f(x,y)=4−14(2x2+y2) at the point (1,2,5/2).

Solution

We define f(x,y) and compute fx(1,2) and fy(1,2).

f [ x _ , y _ ] = 4 − 1 / 4 ( 2 x ∧ 2 + y ∧ 2 ) ;

f [ 1 , 2 ]

dx = D [ f [ x , y ] , x ] /. { x → 1 , y → 2 }

dy = D [ f [ x , y ] , y ] /. { x → 1 , y → 2 }

5 2

−1

−1

Using (3.38), an equation of the tangent plane is z=−1(x−1)−1(y−2)+f(1,2). Using (3.40), parametric equations of the normal line are x=1−t, y=2−t, z=f(1,2)−t. We confirm the result graphically by graphing f(x,y) together with the tangent plane in p1 using Plot3D. We use ParametricPlot3D to graph the normal line in p2 and then display p1 and p2 together with Show in Fig. 3.71.

Figure 3.71 Graph of f(x,y) with a tangent plane and normal line (Washington University in St. Louis colors).

p1 = Plot3D [ f [ x , y ] , { x , − 1 , 3 } , { y , 0 , 4 } ,

PlotStyle → { CMYKColor [ 0 , 1 , . 59 , . 24 ] , Opacity [ . 3 ] } ] ;

p2 = Plot3D [ dx ( x − 1 ) + dy ( y − 2 ) + f [ 1 , 2 ] , { x , − 1 , 3 } , { y , 0 , 4 } ,

PlotStyle → { CMYKColor [ . 59 , . 41 , . 42 , . 14 ] , Opacity [ . 6 ] } ] ;

p3 = ParametricPlot3D [ { 1 + dx t , 2 + dy t , f [ 1 , 2 ] − t } , { t , − 4 , 4 } ,

PlotStyle → { CMYKColor [ 1 , 0 , . 6 , . 4 ] , Thickness [ . 01 ] } ] ;

Show [ p1 , p2 , p3 , PlotRange → { { − 1 , 3 } , { 0 , 4 } , { 0 , 4 } } ,

BoxRatios → Automatic , AxesLabel → { x , y , z } ]

Because z=−1(x−1)−1(y−2)+f(1,2) is the “best” linear approximation of f(x,y) near (1,2), the graphs are very similar near (1,2) as shown in the three-dimensional plot. We also expect the level curves of each near (1,2) to be similar, which is confirmed with ContourPlot in Fig. 3.72.

Figure 3.72 Zooming in near (1,2).

p4 = ContourPlot [ f [ x , y ] , { x , 0.75 , 1.25 } , { y , 1.75 , 2.25 } ,

Contours → 20 ,

ContourStyle → { { CMYKColor [ 0 , 1 , . 59 , . 24 ] , Thickness [ . 01 ] } } ,

PlotLabel → “(a)” ] ;

p5 = ContourPlot [ dx ( x − 1 ) + dy ( y − 2 ) + f [ 1 , 2 ] , { x , 0.75 , 1.25 } , { y , 1.75 , 2.25 } ,

Contours → 20 ,

ContourStyle → { { CMYKColor [ . 59 , . 41 , . 42 , . 14 ] , Thickness [ . 01 ] } } ,

PlotLabel → “(b)” ] ;

Show[GraphicsRow[{p4,p5}]] □

Lagrange Multipliers

Certain types of optimization problems can be solved using the method of Lagrange multipliers that is based on the following theorem.

Theorem 3.19

Lagrange's Theorem

Let z=f(x,y) and z=g(x,y) be real-valued functions with continuous partial derivatives and let z=f(x,y) have an extreme value at a point (x0,y0) on the smooth constraint curve g(x,y)=0. If ▽g(x0,y0)≠0, then there is a real number λ satisfying

▽ f ( x 0 , y 0 ) = λ ▽ g ( x 0 , y 0 ) .

(3.41)

Graphically, the points (x0,y0) at which the extreme values occur correspond to the points where the level curves of z=f(x,y) are tangent to the graph of g(x,y)=0.

Example 3.66

Find the maximum and minimum values of f(x,y)=xy subject to the constraint 14x2+19y2=1.

Solution

For this problem, f(x,y)=xy and g(x,y)=14x2+19y2−1. Observe that parametric equations for 14x2+19y2=1 are x=2cos⁡t, y=3sin⁡t, 0⩽t⩽2π. In Fig. 3.73 (a), we use ParametricPlot3D to parametrically graph g(x,y)=0 and f(x,y) for x and y-values on the curve g(x,y)=0 by graphing

{ x = 2 cos ⁡ t y = 3 sin ⁡ t z = 0 and { x = 2 cos ⁡ t y = 3 sin ⁡ t z = x ⋅ y = 6 cos ⁡ t sin ⁡ t

for 0⩽t⩽2π. Our goal is to find the minimum and maximum values in Fig. 3.73 (a) and the points at which they occur.

Figure 3.73 (a) f(x,y) on g(x,y)=0. (b) Level curves of f(x,y) together with g(x,y)=0 (Yale University colors).

f [ x _ , y _ ] = x y ;

g [ x _ , y _ ] = x ∧ 2 / 4 + y ∧ 2 / 9 − 1 ;

s1 = ParametricPlot3D [ { 2 Cos [ t ] , 3 Sin [ t ] , 0 } , { t , 0 , 2 Pi } ,

PlotStyle → { { CMYKColor [ 1 , . 75 , . 08 , . 4 ] , Thickness [ . 02 ] } } ] ;

s2 = ParametricPlot3D [ { 2 Cos [ t ] , 3 Sin [ t ] , 6 Cos [ t ] Sin [ t ] } , { t , 0 , 2 Pi } ,

PlotStyle → { { CMYKColor [ . 42 , . 4 , . 44 , . 04 ] , Thickness [ . 02 ] } } ] ;

plot1 = Show [ s1 , s2 , BoxRatios → Automatic , PlotRange → All ,

PlotLabel → “(a)” , AxesLabel → { x , y , z } ] ;

To implement the method of Lagrange multipliers, we compute fx(x,y), fy(x,y), gx(x,y), and gy(x,y) with D.

fx = D [ f [ x , y ] , x ]

fy = D [ f [ x , y ] , y ]

gx = D [ g [ x , y ] , x ]

gy=D[g[x,y],y]

y

x

x 2

2 y 9

Solve is used to solve the system of equations (3.41):

f x ( x , y ) = λ g x ( x , y ) f y ( x , y ) = λ g y ( x , y ) g ( x , y ) = 0

for x, y, and λ.

vals = Solve [ { fx = = λ gx , fy = = λ gy , g [ x , y ] = = 0 } , { x , y , λ } ]

{ { x → − 2 , y → − 3 2 , λ → 3 } , { x → − 2 , y → 3 2 , λ → − 3 } , { x → 2 , y → − 3 2 , λ → − 3 } , { x → 2 , y → 3 2 , λ → 3 } }

The corresponding values of f(x,y) are found using ReplaceAll (/.).

n1 = { x , y , f [ x , y ] } /. vals // TableForm

− 2 − 3 2 3 − 2 3 2 − 3 2 − 3 2 − 3 2 3 2 3

N [ n1 ]

− 1.41421 − 2.12132 3 . − 1.41421 2.12132 − 3 . 1.41421 − 2.12132 − 3 . 1.41421 2.12132 3 .

We conclude that the maximum value f(x,y) subject to the constraint g(x,y)=0 is 3 and occurs at (2,322) and (−2,−322). The minimum value is −3 and occurs at (−2,322) and (2,−322). We graph several level curves of f(x,y) and the graph of g(x,y)=0 with ContourPlot and show the graphs together with Show. The minimum and maximum values of f(x,y) subject to the constraint g(x,y)=0 occur at the points where the level curves of f(x,y) are tangent to the graph of g(x,y)=0 as illustrated in Fig. 3.73 (b).

gradfplot = StreamPlot [ Evaluate [ { D [ f [ x , y ] , x ] , D [ f [ x , y ] , y ] } ] , { x , − 3 , 3 } ,

{ y , − 3 , 3 } , StreamStyle → CMYKColor [ 1 , . 75 , . 08 , . 4 ] ] ;

cp1 = ContourPlot [ f [ x , y ] , { x , − 3 , 3 } , { y , − 3 , 3 } , Contours - > 30 ,

ContourShading - > False , PlotPoints - > 40 ,

ContourStyle → { { CMYKColor [ . 42 , . 4 , . 44 , . 04 ] , Thickness [ . 01 ] } } ] ;

cp2 = ContourPlot [ g [ x , y ] == 0 , { x , − 3 , 3 } , { y , − 3 , 3 } , ContourStyle - > Thickness [ 0.01 ] ,

ContourShading - > False ] ;

plot2 = Show [ cp1 , cp2 , gradfplot , Frame → False ,

Axes → Automatic , AxesOrigin → { 0 , 0 } , AxesLabel → { x , y } ,

PlotLabel → “(b)” ] ;

Show [ GraphicsRow [ { plot1 , plot2 } ] ]

Observe that the maximum and minimum values occur where the gradient vectors of z=f(x,y) are parallel to the gradient vectors of z=g(x,y) on the equation g(x,y)=0. □

3.5.3 Iterated Integrals

The Integrate command, used to compute single integrals, is used to compute iterated integrals. The command

Integrate[f[x,y],{y,c,d},{x,a,b}]

attempts to compute the iterated integral

∫ c d ∫ a b f ( x , y ) d x d y .

(3.42)

If Mathematica cannot compute the exact value of the integral, it is returned unevaluated, in which case numerical results may be more useful. The iterated integral (3.42) is numerically evaluated with the command N or

NIntegrate[f[x,y],{y,c,d},{x,a,b}]

Example 3.67

Evaluate each integral: (a) ∫24∫12(2xy2+3x2y)dxdy; (b) ∫02∫y22y(3x2+y3)dxdy; (c) ∫0∞∫0∞xye−x2−y2dydx; (d) ∫0π∫0πesin⁡xydxdy.

Solution

(a) First, we compute ∬(2xy2+3x2y)dxdy with Integrate. Second, we compute ∫24∫12(2xy2+3x2y)dxdy with Integrate.

Integrate [ 2 x y ∧ 2 + 3 x ∧ 2 y , y , x ]

16x2y2(3x+2y)

Integrate [ 2 x y ∧ 2 + 3 x ∧ 2 y , { y , 2 , 4 } , { x , 1 , 2 } ]

98

(b) We illustrate the same commands as in (a), except we are integrating over a nonrectangular region.

Integrate [ 3 x ∧ 2 + y ∧ 3 , { x , y ∧ 2 , 2 y } ]

8 y 3 + 2 y 4 − y 5 − y 6

Integrate [ 3 x ∧ 2 + y ∧ 3 , y , { x , y ∧ 2 , 2 y } ]

2 y 4 + 2 y 5 5 − y 6 6 − y 7 7

Integrate [ 3 x ∧ 2 + y ∧ 3 , { y , 0 , 2 } , { x , y ∧ 2 , 2 y } ]

1664 105

(c) Improper integrals can be handled in the same way as proper integrals.

Integrate [ x   y Exp [ − x ∧ 2 − y ∧ 2 ] , x , y ]

1 4 e − x 2 − y 2

Integrate [ x y Exp [ − x ∧ 2 − y ∧ 2 ] , { x , 0 , Infinity } , { y , 0 , Infinity } ]

1 4

(d) In this case, Mathematica cannot evaluate the integral exactly so we use NIntegrate to obtain an approximation.

Integrate [ Exp [ Sin [ x y ] ] , y , x ]

∫ ∫ e Sin [ x y ] d x d y

NIntegrate [ Exp [ Sin [ x y ] ] , { y , 0 , Pi } , { x , 0 , Pi } ]

15.5092 □

Area, Volume, and Surface Area

Typical applications of iterated integrals include determining the area of a planar region, the volume of a region in three-dimensional space, or the surface area of a region in three-dimensional space. The area of the planar region R is given by

A = ∬ R d A .

(3.43)

If z=f(x,y) has continuous partial derivatives on a closed region R, then the surface area of the portion of the surface that projects onto R is given by

S A = ∬ R ( ∂ f ∂ x ) 2 + ( ∂ f ∂ y ) 2 + 1 d A .

(3.44)

If f(x,y)⩾g(x,y) on R, the volume of the region between the graphs of f(x,y) and g(x,y) is

V = ∬ R ( f ( x , y ) − g ( x , y ) ) d A .

(3.45)

Example 3.68

Find the area of the region R bounded by the graphs of y=2x2 and y=1+x2.

Solution

We begin by graphing y=2x2 and y=1+x2 with Plot in Fig. 3.74. The x-coordinates of the intersection points are found with Solve.

Figure 3.74 y = 2x² and y = 1 + x² for −3/2 ⩽ x ⩽ 3/2 (Michigan State University colors).

Plot [ Tooltip [ { 2 x ∧ 2 , 1 + x ∧ 2 } ] , { x , − 3 / 2 , 3 / 2 } ,

PlotStyle → { { CMYKColor [ . 82 , 0 , . 64 , . 7 ] , Thickness [ . 01 ] } ,

{ CMYKColor [ . 43 , . 3 , . 33 , 0 ] , Thickness [ . 01 ] } } ,

AxesLabel → { x , y } ]

Solve [ 2 x ∧ 2 = = 1 + x ∧ 2 ]

{ { x → − 1 } , { x → 1 } }

Using (3.43) and taking advantage of symmetry, the area of R is given by

A = ∬ R d A = 2 ∫ 0 1 ∫ 2 x 2 1 + x 2 d y d x ,

which we compute with Integrate.

2Integrate[1,{x,0,1},{y,2x∧2,1+x∧2}]

4 3

We conclude that the area of R is 4/3. □

If the problem exhibits “circular symmetry,” changing to polar coordinates is often useful. If R={(r,θ)|a⩽r⩽b,α⩽θ⩽β}, then

∬ R f ( x , y ) d A = ∫ α β ∫ a b f ( r cos ⁡ θ , r sin ⁡ θ ) r d r d θ .

Example 3.69

Find the surface area of the portion of

f ( x , y ) = 4 − x 2 − y 2

that lies above the region R={(x,y)|x2+y2⩽1}.

Solution

First, observe that the domain of f(x,y) is

{ ( x , y ) | − 4 − y 2 ⩽ x ⩽ 4 − y 2 , − 2 ⩽ y ⩽ 2 } = { ( r , θ ) | 0 ⩽ r ⩽ 2 , 0 ⩽ θ ⩽ 2 π } .

Similarly,

R = { ( x , y ) | − 1 − y 2 ⩽ x ⩽ 1 − y 2 , − 1 ⩽ y ⩽ 1 } = { ( r , θ ) | 0 ⩽ r ⩽ 1 , 0 ⩽ θ ⩽ 2 π } .

With this observation, we use ParametricPlot3D to graph f(x,y) in p1 and the portion of the graph of f(x,y) above R in p2 and show the two graphs together with Show. We wish to find the area of the black region in Fig. 3.75.

Figure 3.75 The portion of the graph of f(x,y) above R (Iowa State University colors).

f [ x _ , y _ ] = Sqrt [ 4 − x ∧ 2 − y ∧ 2 ] ;

p1 = ParametricPlot3D [ { r Cos [ t ] , r Sin [ t ] , f [ r Cos [ t ] , r Sin [ t ] ] } , { r , 0 , 2 } ,

{ t , 0 , 2 Pi } , PlotPoints → 45 , PlotStyle → { { CMYKColor [ . 02 , 1 , . 85 , . 06 ] ,

Opacity [ . 4 ] } } ] ;

p2 = ParametricPlot3D [ { r Cos [ t ] , r Sin [ t ] , f [ r Cos [ t ] , r Sin [ t ] ] } , { r , 0 , 1 } ,

{ t , 0 , 2 Pi } , PlotPoints → 45 , PlotStyle → { { CMYKColor [ . 31 , . 38 , . 75 , . 76 ] ,

Opacity [ . 4 ] } } ] ;

Show [ p1 , p2 , BoxRatios → Automatic ]

We compute fx(x,y), fy(x,y) and [fx(x,y)]2+[fy(x,y)]2+1 with D and Simplify.

fx = D [ f [ x , y ] , x ]

fy = D [ f [ x , y ] , y ]

− x 4 − x 2 − y 2

− y 4 − x 2 − y 2

s1 = Simplify [ Sqrt [ 1 + fx ∧ 2 + fy ∧ 2 ] ]

2 − 1 − 4 + x 2 + y 2

Then, using (3.44), the surface area is given by

S A = ∬ R ( ∂ f ∂ x ) 2 + ( ∂ f ∂ y ) 2 + 1 d A = ∬ R 2 4 − x 2 − y 2 d A = ∫ − 1 1 ∫ − 1 − y 2 1 − y 2 2 4 − x 2 − y 2 d x d y .

(3.46)

However, notice that in polar coordinates,

R = { ( r , θ ) | 0 ⩽ r ⩽ 1 , 0 ⩽ θ ⩽ 2 π }

so in polar coordinates the surface area is given by

S A = ∫ 0 2 π ∫ 0 1 2 4 − r 2 r d r d θ ,

s2 = Simplify [ s1 /. { x → r Cos [ t ] , y → r Sin [ t ] } ]

2 1 4 − r 2

which is much easier to evaluate than (3.46). We evaluate the iterated integral with Integrate

s3=Integrate[rs2,{t,0,2Pi},{r,0,1}]

− 4 ( − 2 + 3 ) π

N [ s3 ]

3.36715

and conclude that the surface area is (8−43)π≈3.367. □

Example 3.70

Find the volume of the region between the graphs of z=4−x2−y2 and z=2−x.

Solution

We begin by graphing z=4−x2−y2 and z=2−x together with Plot3D in Fig. 3.76 (a). The region of integration, R, is determined by graphing 4−x2−y2=2−x with ContourPlot in Fig. 3.76 (b).

Figure 3.76 (a) z = 4 − x² − y² and z = 2 − x for −2 ⩽ x ⩽ 2 and −2 ⩽ y ⩽ 2. (b) Graph of 4 − x² − y² = 2 − x (University of California-Irvine colors).

p1 = Plot3D [ { 4 − x ∧ 2 − y ∧ 2 , 2 − x } , { x , − 2 , 2 } , { y , − 2 , 2 } ,

PlotRange → { { − 2 , 2 } , { − 2 , 2 } , { − 2 , 4 } } ,

PlotStyle → { { CMYKColor [ . 93 , . 62 , . 09 , . 01 ] , Opacity [ . 4 ] } ,

{ { CMYKColor [ . 01 , . 16 , 1 , 0 ] , Opacity [ . 6 ] } } } ,

MeshFunctions - > { # 3 & } ,

BoxRatios → Automatic ] ;

p2 = ContourPlot [ 4 − x ∧ 2 − y ∧ 2 − ( 2 − x ) = = 0 , { x , − 2 , 2 } , { y , − 2 , 2 } ,

PlotPoints → 50 , Frame → False , Axes → Automatic , AxesOrigin → { 0 , 0 } ,

ContourStyle → { { CMYKColor [ . 01 , . 16 , 1 , . 0 ] , Thickness [ . 01 ] } } ] ;

Show[GraphicsRow[{p1,p2}]]

Another way to see the situation illustrated in Fig. 3.77 is to use RegionPlot3D, which works in the same way as RegionPlot but in three dimensions.

Figure 3.77 Using RegionPlot3D to see the region 2 − x ⩽ z ⩽ 4 − x² − y² (University of California-Irvine colors).

RegionPlot3D [ 2 − x ⩽ z < = 4 − x ∧ 2 − y ∧ 2 , { x , − 2 , 2 } , { y , − 2 , 2 } ,

{ z , − 4 , 4 } , PlotPoints → 75 , Mesh → False ,

PlotStyle - > { CMYKColor [ . 93 , . 62 , . 09 , . 01 ] , Opacity [ . 4 ] } ]

Completing the square shows us that

R = { ( x , y ) | ( x − 1 2 ) 2 + y 2 ⩽ 9 4 } = { ( x , y ) | 1 2 − 1 2 9 − 4 y 2 ⩽ x ⩽ 1 2 + 1 2 9 − 4 y 2 , − 3 2 ⩽ y ⩽ 3 2 } .

Thus, using (3.45), the volume of the solid is given by

V = ∬ R [ ( 4 − x 2 − y 2 ) − ( 2 − x ) ] d A = ∫ − 3 2 3 2 ∫ 1 2 − 1 2 9 − 4 y 2 1 2 + 1 2 9 − 4 y 2 [ ( 4 − x 2 − y 2 ) − ( 2 − x ) ] d x d y ,

which we evaluate with Integrate.

i1 = Integrate [ ( 4 − x ∧ 2 − y ∧ 2 ) − ( 2 − x ) , { y , − 3 / 2 , 3 / 2 } ,

{x,1/2−1/2Sqrt[9−4y∧2],1/2+1/2Sqrt[9−4y∧2]}]

81 π 32

N [ i1 ]

7.95216

We conclude that the volume is 8132π≈7.952. □

Triple Iterated Integrals

Triple iterated integrals are calculated in the same manner as double iterated integrals.

Example 3.71

Evaluate

∫ 0 π / 4 ∫ 0 y ∫ 0 y + z ( x + 2 z ) sin ⁡ y d x d z d y .

Solution

Entering

i1 = Integrate [ ( x + 2 z ) Sin [ y ] , { y , 0 , Pi / 4 } , { z , 0 , y } , { x , 0 , y + z } ]

− 17 ( 384 − 96 π − 12 π 2 + π 3 ) 384 2

calculates the triple integral exactly with Integrate.

An approximation of the exact value is found with N.

N [ i1 ]

0.157206 □

We illustrate how triple integrals can be used to find the volume of a solid when using spherical coordinates.

Example 3.72

Find the volume of the torus with equation in spherical coordinates ρ=sin⁡ϕ.

Solution

We proceed by graphing the torus with SphericalPlot3D in Fig. 3.79 (see Fig. 3.78 for the help feature associated with this command).

Figure 3.78 Mathematica's help for SphericalPlot3D.

Figure 3.79 A graph of the torus (University of Arizona colors).

SphericalPlot3D [ Sin [ Phi ] , { Phi , 0 , Pi } , { theta , 0 , 2 Pi } , PlotPoints → 40 ]

In general, the volume of the solid region D is given by

V = ∭ D d V .

Thus, the volume of the torus is given by the triple iterated integral

V = ∫ 0 2 π ∫ 0 π ∫ 0 sin ⁡ ϕ ρ 2 sin ⁡ ϕ d ρ d ϕ d θ ,

i1 = Integrate [ rho ∧ 2 Sin [ phi ] , { theta , 0 , 2 Pi } ,

{ phi , 0 , Pi } , { rho , 0 , Sin [ phi ] } ]

π24

N [ i1 ]

2.4674

which we evaluate with Integrate. We conclude that the volume of the torus is 14π2≈2.467. □