MATLAB Differential and Integral Calculus
MATLAB Differential and Integral Calculus
MATLAB Differential and Integral Calculus
com
Contents at a Glance
iii
www.pdfgrip.com
Chapter 1
1
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
4. We can also work with complex numbers. We find the result of the operation raising
(2 + 3i) to the power 10 by typing the expression (2 + 3i) ^ 10.
>> (2 + 3i) ^ 10
ans =
-1 415249999999998e + 005 - 1. 456680000000000e + 005i
5. The previous result is also available in short format, using the “format short” command.
>> format short
>> (2 + 3i)^10
ans =
-1.4152e + 005- 1.4567e + 005i
6. We can calculate the value of the Bessel function J0 at 11.5. To do this we type
besselj(0,11.5).
>> besselj(0,11.5)
ans =
-0.0677
7. We can also perform numerical integration. To calculate the integral of sin(sin(x)) between
0 and p we type int(‘sin((sin(x))’, 0, pi).
>> int ('sin(sin(x))', 0, pi)
ans =
1235191162052677/2251799813685248 * pi
These ideas will be treated more thoroughly later in the book.
2
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
1. We can expand the following algebraic expression: ((x + 1)(x + 2) - (x + 2) ^ 2)^3. This
is done by typing: expand(‘((x + 1)(x + 2) - (x + 2) ^ 2) ^ 3’). The result will be another
algebraic expression:
>> syms x; expand(((x + 1) *(x + 2)-(x + 2) ^ 2) ^ 3)
ans =
-x ^ 3-6 * x ^ 2-12 * x-8
2. We can factor the result of the calculation in the above example by typing:
factor(‘((x + 1) *(x + 2) - (x + 2) ^ 2) ^ 3’)
>> syms x; factor(((x + 1)*(x + 2)-(x + 2)^2)^3)
ans =
-(x+2)^3
3. We can find the indefinite integral of the function (x ^ 2) sin(x) ^ 2 by typing:
int(‘x ^ 2 * sin(x) ^ 2’, ‘x’)
>> int('x^2*sin(x)^2', 'x')
ans =
x ^ 2 *(-1/2 * cos(x) * sin(x) + 1/2 * x)-1/2 * x * cos(x) ^ 2 + 1/4 * cos(x) * sin(x) +
1/4 * 1/x-3 * x ^ 3
4. We can simplify the previous result:
>> syms x; simplify(int(x^2*sin(x)^2, x))
ans =
sin(2*x)/8 -(x*cos(2*x))/4 -(x^2*sin(2*x))/4 + x^3/6
5. We can present the previous result using a more elegant mathematical notation:
>> syms x; pretty(simplify(int(x^2*sin(x)^2, x)))
ans =
2 3
sin(2 x) x cos(2 x) x sin(2 x) x
-------- - ---------- - ----------- + --
8 4 4 6
3
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
6. We can find the series expansion up to order 12 of the function x ^ 2 * sin (x) ^ 2,
presenting the result in elegant form:
>> pretty(taylor('x^2*sin(x)^2',12))
4 6 8 10 12
x - 1/3 x + 2/45 x - 1/315 x + o (x)
7. We can solve the equation 3ax - 7 x ^ 2 + x ^ 3 = 0 (where a is a parameter):
>> solve('3*a*x-7*x^2 + x^3 = 0', 'x')
ans =
[ 0]
[7/2 + 1/2 *(49-12*a) ^(1/2)]
[7/2-1/2 *(49-12*a) ^(1/2)]
8. We can find the five solutions of the equation x ^ 5 + 2 x + 1 = 0:
>> solve('x^5+2*x+1','x')
ans =
RootOf(_Z^5+2*_Z+1)
As the result does not explicitly give five solutions, we apply the “allvalues” command:
>> allvalues(solve('x^5+2*x+1','x'))
ans =
[-.7018735688558619-. 8796971979298240 * i]
[-. 7018735688558619 +. 8796971979298240 * i]
[-. 4863890359345430]
[.9450680868231334-. 8545175144390459 * i]
[. 9450680868231334 +. 8545175144390459 * i]
On the other hand, MATLAB can use the Maple program libraries to work with symbolic math, and can thus
extend its field of action. In this way, MATLAB can be used to work on such topics as differential forms, Euclidean
geometry, projective geometry, statistics, etc.
At the same time, Maple can also benefit from MATLAB’s powers of numerical calculation, which might be used,
for example, in combination with the Maple libraries (combinatorics, optimization, number theory, etc.)
4
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-1.
2. We can give the above graph a title and label the axes, and we can add a grid. See Figure 1-2.
>> x = linspace(-pi/4,pi/4,300);
>> y=x.*sin(1./x);
>> plot(x,y);
>> grid;
>> xlabel('Independent variable X');
>> ylabel('Dependent variable Y');
>> title('The function y=xsin(1/x)')
6
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-2.
7
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-3.
These 3D graphics allow you to get a clear picture of figures in space, and are very helpful in visually identifying
intersections between different bodies, and in generating all kinds of space curves, surfaces and volumes of revolution.
4. We can generate the three dimensional graph corresponding to the helix with parametric
coordinates: x = sin(t), y = cos(t), z = t. See Figure 1-4.
>> t=0:pi/50:10*pi;
>> plot3(sin(t),cos(t),t)
Figure 1-4.
8
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
5. We can represent a planar curve given by its polar coordinates r = cos(2t) * sin(2t) for
t varying in the range between 0 and p by equally spaced points 0.01 apart. See Figure 1-5.
>> t = 0:. 1:2 * pi;
>> r = sin(2*t). * cos(2*t);
>> polar(t,r)
Figure 1-5.
6. We can make a graph of a symbolic function using the command “ezplot”. See Figure 1-6.
>> y ='x ^ 3 /(x^2-1)';
>> ezplot(y,[-5,5])
Figure 1-6.
9
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-7.
At other times, depending on the type of entry (user input) given to MATLAB, the response is returned using the
expression “ans =”. See Figure 1-8.
10
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-8.
It is important to pay attention to the use of uppercase versus lowercase letters, parentheses versus square
brackets, spaces and punctuation (particularly commas and semicolons).
11
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
12
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
14
www.pdfgrip.com
Chapter 1 ■ Introduction and the MATLAB Environment
Figure 1-9.
The command >>dos_command is used to execute the DOS command in the MATLAB screen. Using the three
previous commands, not only DOS commands, but also all kinds of executable files or batch tasks can be executed.
The command >>dos dos_command is also used to execute the specified DOS command in automatic mode in
the MATLAB Command Window.
To exit MATLAB, simply type quit in the Command Window, and then press Enter.
Chapter 2
MATLAB provides commands that allow you to calculate virtually all types of limits. The same functions are used to
calculate limits of sequences and limits of functions. The commands for the analysis of one and several variables are
similar. In this chapter we will present numerous exercises which illustrate MATLAB’s capabilities in this field.
The syntax of the commands concerning limits are presented below:
maple (‘limit(sequence, n=infinity)’) or limit(sequence, n, inf) or limit(sequence, inf)
calculates the limit as n tends to infinity of the sequence defined by its general term.
maple (‘limit(function, x=a)’) or limit(function, x, a) or limit(function, a) calculates
the limit of the function of the variable x, indicated by its analytical expression, as the
variable x tends towards the value a.
maple (‘limit(function, x=a, right)’) or limit(function, x, a, ‘right’) calculates the limit
of the function of the variable x, indicated by its analytical expression, as the variable x
tends to the value a from the right.
maple (‘limit(function, x=a, left)’) or limit(function, x, a, ‘left’) calculates the limit of
the function of the variable x, indicated by its analytical expression, as the variable x
tends to the value a from the left.
maple (‘limit(expr, var=a, complex)’) is the complex limit of expr as the variable var
tends to the value a.
maple (‘limit(expr, {v1=a1, …, vn=an})’) is the n-dimensional limit of expr as v1 tends
to a1, v2 tends to a2,…, vn tends to an.
maple (‘Limit(expr, var=a)’) or maple (‘Limit(expr, {v1=a1, …, vn=an})’) is the inert
limit of the expression expr for the specified values of the variable or variables.
17
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-1
æ -3 + 2n ö
4
1 + 7n 2 + 3n 3 é 1 + n ö4 1 + n ù 1+ n
lim ç ÷ , lim , lim êæç ÷ 5 ú
, lim n 2
n®¥ -7 + 3n n®¥ 5 - 8n + 4n 3
è ø n®¥
ëêè 2 ø n úû n®¥ n
In the first two limits we face the typical uncertainty given by the quotient ¥ ¥ :
>> syms n
>> limit(((2*n-3)/(3*n-7))^4, inf)
ans =
16/81
>> limit((3*n^3+7*n^2+1)/(4*n^3-8*n+5),n,inf)
ans =
3/4
The last two limits present an uncertainty of the form ∞.0 and ∞0:
>> limit(((n+1)/2)*((n^4+1)/n^5), inf)
ans =
1/2
>> limit(((n+1)/n^2)^(1/n), inf)
ans =
1
18
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-2
19
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-3
sin é( ax ) ù
2
x
-1 - x x - 2+x x ë û , lim e - 1 .
lim , lim , lim 1 + x , lim
x ®1 -1 + x x ®2 -3 + 1 + 4 x x ®0 x ®0 x 2 x ®0 log (1 + x )
Initially, we have two indeterminates of type 0/0 and one of the form 1∞:
>> syms x
>> limit((x-1)/(x^(1/2)-1),x,1)
2
>> limit((x-(x+2)^(1/2))/((4*x+1)^(1/2)-3),2)
9/8
>> limit((1+x)^(1/x))
exp (1)
The last two are indeterminates of the form 0/0:
>> syms x a, limit(sin(a*x)^2/x^2,x,0)
a^2
>> numeric(limit((exp(1)^x-1)/log(1+x)))
1
20
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-4
Figure 2-1.
21
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
By simple observation, we see that the limit from the right is 1 and the limit from the left is - 1.
For the next two limits we have:
>> limit(abs(x^2-x-7),x,3)
ans =
1
>> limit((x-1)/(x^n-1),x,1)
ans =
1/n
For the last limit we have the following:
>> limit(exp(1)^(1/x),x,0)
ans =
NaN
>> limit(exp(1)^(1/x),x,0,'left')
ans =
0
>> limit(exp(1)^(1/x),x,0,'right')
ans =
INF
Then, there is no limit at x = 0. If we plot the function (see Figure 2-2), we have:
>> ezplot (exp(1)^(1/x), [- 10, 10 - 3, 3])
22
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-2.
We see that the function becomes (positive) infinite at 0 when approaching from the right, and it tends to 0 when
approaching from the left. Thus we conclude that the function has no limit at x = 0.
Exercise 2-5
23
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
>> pretty(sym(maple('limit(x^n/(1+x^n),n=infinity)')))
0
>> pretty(sym(maple('limit(x^n/(n+x^n),n=infinity) ')))
0
>> pretty(sym(maple('limit(piecewise(x>=1/(n+1) and x<=1/n,
sin(Pi/x)^2,x<1/(n+1) and x<1/n,0),n=infinity) ')))
0
Now we support these results with the graphical representations shown in Figures 2-3, 2-4 and 2-5, which in
each case illustrates the convergence of the sequence of functions towards the constant 0 function. The
syntax for the plots is as follows:
>> fplot('[x,x^2/2,x^3/3,x^4/4,x^5/5,x^6/6,x^7/7,x^8/8,x^9/9,x^10/10]',[0,1,-1/2,1])
>> fplot('[x/(1+x),x^2/(1+x^2),x^3/(1+x^3),x^4/(1+x^4),x^5/(1+x^5),x^6/(1+x^6),
x^7/(1+x^7),x^8/(1+x^8),x^9/(1+x^9),x^10/(1+x^10)]',[0,1,-1/2,1])
>> fplot('[x/(1+x),x^2/(2+x^2),x^3/(3+x^3),x^4/(4+x^4),x^5/(5+x^5),x^6/(6+x^6),
x^7/(7+x^7),x^8/(8+x^8),x^9/(9+x^9),x^10/(10+x^10)]',[0,1,-1/2,1])
Figure 2-3.
24
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-4.
Figure 2-5.
Exercise 2-6
25
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
26
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-6.
Figure 2-7.
27
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-8.
We keep the definitions of the sequences of functions in the M-files named seq1.m, seq2.m and seq3.m
respectively. The graphical respresentation of the first 10 functions of each sequence is as follows:
>> fplot ('[seq1(x,1), seq1(x,2), seq1(x,3), seq1(x,4), seq1(x,5),])
([seq1(x,6), seq1(x,7), seq1(x,8), seq1(x,9), seq1(x,10)]', [- 2, 2, - 2, 2])
>> fplot ('[seq2(x,1), seq2(x,2), seq2(x,3), seq2(x,4), seq2(x,5),])
([seq2(x,6), seq2(x,7), seq2(x,8), seq2(x,9), seq2(x,10)]'[- 10, 10, 0, 3/2])
>> fplot ('[seq3(x,1), seq3(x,2), seq3(x,3), seq3(x,4), seq3(x,5),])
([seq3(x,6), seq3(x,7), seq3(x,8), seq3(x,9), seq3(x,10)]', [- 0, 1, - 1/2, 10])
2.4 Continuity
A function f is continuous at the point x = a if:
lim f ( x ) = f ( a ) .
x ®a
Otherwise, it is discontinuous at the point. In other words, in order for a function f to be continuous at a it must
be defined at a and the limit of the function at a must be equal to the value of the function at a.
If the limit of f (x) as x tends to a exists but is different to f (a), then f is discontinuous at a, and we say f has an
avoidable discontinuity at a. The discontinuity is resolved by redefining f (a) to coincide with the limit.
If the two lateral limits of f (x) at a exist (whether finite or infinite) but are different, then the discontinuity of f at a
is said to be of the first kind. The difference between the two lateral limits is the jump. If the jump is finite then the
discontinuity is said to be of the first kind with finite jump, otherwise it is of the first kind with infinite jump.
If either of the lateral limits do not exist, then the discontinuity is said to be of the second kind.
28
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-7
The problem arises at the point x = 0, at which the function f is not defined. Therefore, the function is
discontinuous at x = 0, This discontinuity can be avoided by redefining the function at x = 0 with a value
equal to lim f ( x ) .
x ®0
>> limit(sin(x)/x,x,0)
ans =
1
Thus we conclude that the function f ( x ) = sin( x ) / x presents an avoidable discontinuity at x = 0 that is avoided
by defining f (0) = 1. The function is continuous at all non-zero points.
The function g is continuous at any non-zero point a, so lim g ( x ) = g ( a ) .
x ®a
>> limit (sin (1/x), x, a)
ans =
sin(1/a)
29
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
The problem arises at the point x = 0, where the function g is not defined. Therefore, the function is discontinuous
at x = 0. To try to avoid the discontinuity, we calculate lim g ( x ) .
x ®0
>> limit(sin(1/x),x,0)
ans =
-1 .. 1
>> limit(sin(1/x),x,0,'left')
ans =
-1 .. 1
>> limit(sin(1/x),x,0,'right')
ans =
-1 .. 1
We see that the limit does not exist at x = 0 (the limit has to be unique and here the result given is all points in
the interval [− 1,1]), and neither of the lateral limits exist. Thus the function has a discontinuity of the second
kind at x = 0.
MATLAB responded to the calculation of the above limits with the expression "– 1.. 1". This is because the graph
of g (x) presents infinitely many oscillations between – 1 and 1. This is illustrated in Figure 2-9:
>> fplot ('sin (1/x)', [- 1, 1])
Figure 2-9.
30
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-8
The only problematic point is x = 0. Now, the function does exist at x = 0 (it has the value 1). We will try to find
the lateral limits as x tends to 0:
>> syms x
>> limit(1/(1+exp(1/x)),x,0,'right')
ans =
0
>> limit(1/(1+exp(1/x)),x,0,'left')
ans =
1
As the lateral boundaries are different, the limit of the function at 0 does not exist. However, the lateral boundaries
are both finite, so the discontinuity at 0 is of the first kind with a finite jump. We illustrate this result in Figure 2-10.
>> fplot('1/(1+exp(1/x))',[-5,5])
Figure 2-10.
31
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-9
The only problematic point is x = 0. The function is defined at x = 0 (it has the value 1). We will try to find the
lateral limits at 0:
>> limit ((exp (x)), x, 0, 'right')
ans =
inf
>> limit((exp(1/x)),x,0, 'left')
ans =
0
As the lateral boundaries are different, the limit of the function as x tends to 0 does not exist. As the right
lateral limit is infinite, the discontinuity of the first kind at x = 0 is an infinite jump. We illustrate this result in
Figure 2-11 above.
>> fplot ('exp (x)', [- 150, 150])
Figure 2-11.
32
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
This characterization theorem allows us to calculate limits of sequences of points in m-dimensional space.
There is another theorem, similar to the above, which characterizes the limits of functions between spaces of
more than one dimension. This theorem enables us to calculate the limits of multivariable functions.
If f : Rn → Rm is a function whose m components are (f1, …, fm). Then it follows that:
lim
x1 ®a1 , x 2 ®a2 ,¼,xn ®an
( f ( x ,x
1 1 2 ,¼, x n ) , f 2 ( x1 , x 2 ,¼, x n ) ,¼, fm ( x1 , x 2 ,¼, x n ) )
= (l1 , l2 ,¼, lm )
if and only if
lim
x1 ®a1 , x 2 ®a2 ,¼,xn ®an
( f (x ,x
1 1 2 ,¼, x n ) ) = l1 ,
lim
x1 ®a1 , x 2 ®a2 ,¼,xn ®an
( f (x ,x
2 1 2 ,¼, x n ) ) = l2 ,
…,
lim
x1 ®a1 , x 2 ®a2 ,¼,xn ®an
( f ( x ,x
m 1 2 ,¼, x n ) ) = lm .
Exercise 2-10
é 1 + n æ 1 ö2 n n ù
lim ê , ç1 + ÷ , ú
n®¥
êë n è n ø 2n - 1 úû
>> pretty(sym(maple('vector([limit((n+1)/n,n=infinity),limit((1+1/n)^(2*n),
n=infinity),limit(n/(2*n-1),n=infinity)])')))
[1-exp (2) 1/2]
33
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-11
é n 1 n 1 + n2 ù
lim ê n , n , 5n , 2 ú
n®¥
ë 1 + n2 n n û
>> pretty(sym(maple('vector([limit((n/(n^2+1))^(1/n),n=infinity),limit((1/n)^(1/n),
n=infinity),limit((5*n)^(1/n) ,n=infinity),limit((n^2+1)/n^2,n=infinity)])')))
[1, 1, 1, 1]
Exercise 2-12
æ sin( x ) x ö
f (x) = ç , 1 + x ÷.
è x ø
>> pretty(sym(maple('vector([limit(sin(x)/x,x=0),limit((1+x)^(1/x),x=0)])')))
[1, exp (1)]
Exercise 2-13
34
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
è x2 ®a2
x1 ®a1 (
lim æç lim ¼ lim f ( x1 , x 2 ,¼, x n ) ö÷
xn ®an ø )
or any of the other limits obtained by permuting the order of the component limits.
The directional limit of f at the point (a1, …, an) depends on the direction of the curve h (t) = (h1(t), h2(t), …, hn(t)),
where h(t0) = (a1, a2, …, an), and is defined to be the value:
lim( f (h(t )) = lim f ( x1 , x 2 ,¼, x n )
t ®t 0 ( x1 ,x 2 ,¼,xn )®( a1 ,a2 ,¼,an )
A necessary condition for a function of several variables to have a limit at a point is that all the iterated limits have
the same value (which will be equal to the value of the limit of the function, if it exists).
It can also happen that the directional limit of a function will vary according to the curve used, so that a different
limit exists for different curves, or the limit exists for some curves and not others.
Another necessary condition for a function of several variables to have a limit at a point is that all directional
limits, i.e. the limit for all curves, have the same value.
Therefore, to prove that a function has no limit at a point it is enough to show that either an iterated limit does
not exist or that two iterated limits have a different value, or we can show that a directional limit does not exist or that
two directional limits have different values.
A practical procedure for calculating the limit of a function of several variables is to change from cartesian to
polar coordinates.
Exercise 2-14
xy
f (x , y) =
x2 + y2
>> syms x y
>> limit (limit ((x*y) /(x^2+y^2), x, 0), y, 0)
ans =
0
>> limit (limit ((x*y) /(x^2+y^2), y, 0), x, 0)
ans =
0
35
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Thus the two iterated limits are the same. Next we calculate the directional limits corresponding to the family of
straight lines y = mx:
>> syms m
>> limit((m*x^2)/(x^2+(m^2)*(x^2)),x,0)
ans =
m /(1+m^2)
The directional limits depend on the parameter m, which will be different for different values of m (corresponding
to different straight lines). Thus, we conclude that the function has no limit at (0,0).
Exercise 2-15
( y 2 - x 2 )2
f (x , y) =
x2 + y4
>> syms x y
>> limit (limit ((y^2-x^2) ^ 2 /(y^4+x^2), x, 0), y, 0)
ans =
1
>> limit (limit ((y^2-x^2) ^ 2 /(y^4+x^2), y, 0), x, 0)
ans =
0
As the two iterated limits are different, we conclude that the function has no limit at the point (0,0).
36
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-16
( y 2 - x )2
f (x , y) =
x2 + y4
>> syms x y m, limit (limit ((y^2-x) ^ 2 /(y^4+x^2), y, 0), x, 0)
ans =
1
>> limit (limit ((y^2-x) ^ 2 /(y^4+x^2), x, 0), y, 0)
ans =
1
Thus the two iterated limits are the same. Next we calculate the directional limits corresponding to the family
of straight lines y = mx:
>> limit(((m*x)^2-x)^2/((m*x)^4+x^2),x,0)
ans =
1
The directional limits corresponding to the family of straight lines y = mx do not depend on m and coincide with
the iterated limits. Next we find the directional limits corresponding to the family of parabolas y ^ 2 = mx:
>> limit(((m*x)-x)^2/((m*x)^2+x^2),x,0)
ans =
(m-1) ^ 2-/(m^2+1)
Thus the directional limits corresponding to this family of parabolas depend on the parameter m, so they are
different. This leads us to conclude that the function has no limit at (0,0).
37
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-17
x2y
f (x , y) =
x2 + y2
>> syms x y, limit (limit ((x^2*y) /(x^2+y^2), x, 0), y, 0)
ans =
0
>> limit (limit ((x^2*y) /(x^2+y^2), x, 0), y, 0)
ans =
0
>> limit(((x^2)*(m*x))/(x^2+(m*x)^2),x,0)
ans =
0
>> limit (((m*y) ^ 2) * y / ((m*y) ^ 2 + y ^ 2), y, 0)
ans =
0
We see that the iterated limits and directional limits corresponding to the given family of lines and parabolas
coincide and are all zero. This leads us to suspect that the limit of the function may be zero. To confirm this,
we transform to polar coordinates and find the limit:
>> syms a r, limit (limit (((r^2) * (cos (a) ^ 2) * (r) * (sin (a))) / ((r^2) *
(cos (a) ^ 2)+(r^2) * (sin (a) ^ 2)), r, 0), a, 0)
ans =
0
Therefore we conclude that the limit of the function is zero at the point (0,0).
This is an example where, as a last resort, we had to transform to polar coordinates. In the above examples
we used families of lines and parabolas, but other curves can be used. The change to polar coordinates can
be crucial in determining limits of functions of several variables. As we have seen, there are sufficient criteria
to show that a function has no limit at a point. However, we do not have necessary and sufficient conditions to
ensure the existence of the limit.
38
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-18
( x - 1)2 y 2
f (x , y) =
( x - 1)2 + y 2
>> syms x y m a r
>> limit (limit (y ^ 2 *(x-1) ^ 2 / (y ^ 2 +(x-1) ^ 2), x, 0), y, 0)
ans =
0
>> limit (limit (y ^ 2 *(x-1) ^ 2 / (y ^ 2 +(x-1) ^ 2), y, 0), x, 0)
ans =
0
>> limit((m*x)^2*(x-1)^2/((m*x)^2+(x-1)^2),x,0)
ans =
0
>> limit((m*x)*(x-1)^2/((m*x)+(x-1)^2),x,0)
ans =
0
We see that the iterated and directional limits coincide. We calculate the limit by converting to polar coordinates:
>> limit (limit ((r ^ 2 * sin (a) ^ 2) * (r * cos (a) - 1) ^ 2 / ((r ^ 2 * sin (a) ^ 2) +
(r * cos (a) - 1) ^ 2), r, 1), a, 0)
ans =
0
The limit is zero at the point (1,0). The surface is depicted in Figure 2-12 where we see the function tends to 0 in
a neighborhood of (1,0):
>> [x, y] = meshgrid(0:0.05:2,-2:0.05:2);
>> z=y.^2.*(x-1).^2./(y.^2+(x-1).^2);
>> mesh(x,y,z), view ([- 23, 30])
39
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-12.
The only problematic point is the origin. We will analyze the continuity of the function at the origin by calculating
its limit there.
40
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
>> syms x y m a r
>> limit (limit ((x^3+y^3) /(x^2+y^2), x, 0), y, 0)
ans =
0
>> limit (limit ((x^3+y^3) /(x^2+y^2), y, 0), x, 0)
ans =
0
>> limit((x^3+(m*x)^3)/(x^2+(m*x)^2),x,0)
ans =
0
We see that the iterated and the linear directional limits coincide. We try to calculate the limit by converting to
polar coordinates:
>> maple ('limit (limit (((r * cos (a)) ^ 3 + (r * sin (a)) ^ 3) / ((r * cos (a)) ^ 2 +
(r * sin (a)) ^ 2), r = 0), a = 0)')
ans =
0
We see that the limit at (0,0) coincides with f(0,0), so the function is continuous at (0,0). The graph in Figure 2-13
confirms the continuity.
>> [x, y] = meshgrid(-1:0.007:1,-1:0.007:1);
>> z=(x.^3+y.^3)./(x.^2+y.^2);
>> mesh (x, y, z)
>> view ([50, - 20])
41
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-13.
Exercise 2-20
42
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Exercise 2-21
f (x, y ) =
[1 - cos( x )]sin( y ) if ( x , y ) ¹ (0 , 0) and f ( 0 , 0 ) = 0.
x +y
3 3
43
www.pdfgrip.com
Chapter 2 ■ Limits and Continuity. One and Several Variables
Figure 2-14.
44
www.pdfgrip.com
Chapter 3
This chapter demonstrates the wide range of features that MATLAB offers which can be used to treat numerical series.
These include the determination of the radius of convergence of a power series, summation of convergent series,
alternating series and so on.
Again, if the limit is 1, we cannot say whether the series diverges or converges.
If a limit of 1 is obtained in both the ratio and root test, we can often use the criterion of Raabe or Duhamel,
which reads as follows:
¥
é æ a(n + 1) ö ù
å a(n) is convergent if lim ên ç 1 - ÷ú > 1
a(n) ø û
n =1
n®¥
ë è
45
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
¥
é æ a(n + 1) ö ù
å a(n) is divergent if lim ên ç 1 - a(n) ø û
÷ú < 1
n =1
n®¥
ë è
However, if the limit is 1, we still cannot conclude anything about the convergence or divergence of the series.
Another useful criterion is the following:
Other criteria such as Gauss majorization, comparison tests and so on can also be implemented via Maple. We
will see some examples below:
EXERCISE 3-1
Study the convergence and, if possible, find the sum of the following series:
¥
1+ n
å n(2 + n)(3 + n)
n =1
46
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-2
Study the convergence and, if possible, find the sum of the following series:
¥ ¥ ¥
n nn nn
å2
n =1
n
, å n! , å 3
n =1 n =1
n
n!
47
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-3
Study the convergence and, if possible, find the sum of the following series:
¥
3 + 2n ¥
n ¥
( n + p )! , p = parameter
å7
n =1 n (1 + n)
n
, åp
n=1
n
, å
n =1
n
p n!p!
48
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
49
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-4
Study the convergence and, if possible, find the sum of the following series:
-n
¥
1 2 ¥ éæ 1 + n ö1+n 1 + n ù
å (1 + ) - n , å êç ÷ - ú
n =1 n n =1 êëè n ø n úû
50
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-5
Study the convergence and, if possible, find the sum of the following series:
n2
¥
5 ¥ ¥
æ n 2 + 2n + 1 ö
å , å( n
n - 1) ,
n
å ç 2 ÷
n=1 2n n =1 n =1 è n + n -1 ø
51
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-6
Study the convergence and, if possible, find the sum of the following series:
¥
æ qö
å tan
n =1
n
çp + ÷×
è nø
EXERCISE 3-7
Study the convergence and, if possible, find the sum of the following series:
¥ ¥
1 1
å n Log (n) , å n [ Log (n)]
n =1 n =1
p
, p = parameter > 0
52
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
(- p) (- p)
n log(2)
When p < 1, this series dominates the series with general term n –p = 1/np. Thus the initial series also diverges.
When p < 1, this series is dominated by the convergent series with general term n –p = 1/n p. Thus the initial
series also converges.
When p = 1, the series reduces to the series studied above, i.e. it diverges.
MATLAB does not offer the sum of either of the two series of this problem.
EXERCISE 3-8
Study the convergence and, if possible, find the sum of the following series:
¥
(n + 1)(n + 2) ¥
1
å , å
(1 + n )
2
n =1 n5 n =1
We will begin by studying the second series. We try to apply the ratio, Raabe and root tests:
>> maple('a:=n->1/(1+n^(1/2))^2');
>> maple('limit(a(n+1)/a(n), n=infinity)')
ans =
1
>> maple('limit(n*(1-a(n+1)/a(n)),n=infinity)')
ans =
1
>> maple('limit(simplify((a(n))^(1/n)), n=infinity)')
ans =
1
Thus all the limits are 1. Therefore, at the moment, we cannot conclude anything about the convergence of the series.
We now compare our series with the divergent harmonic series by finding the limit of the quotient of the
respective general terms:
>> maple('limit(a(n)/(1/n),n=infinity)')
ans =
1
53
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
As the limit is greater than zero, the initial series is also divergent.
We will now analyze the first series of the problem directly, comparing it with the convergent series with general
term 1/n3 by examining the limit of the quotient of the general terms:
>> maple('a:=n->(n+1)*(n+2)/n^5');
>> maple('limit(a(n)/(1/(n^3)),n=infinity)')
ans =
1
As the limit is greater than 0, the initial series is also convergent:
>> maple('sum(a(n),n=1..+infinity)')
ans =
2 * zeta (5) + 1/30 * pi ^ 4 + zeta (3)
Now, we try to approximate the result:
>> maple('evalf(sum(a(n),n=1..+infinity))')
ans =
6.522882114579747
We could have tried to determine whether the first series converges or diverges by using the ratio, root and
Raabe criteria:
>> maple('limit(a(n+1)/a(n), n=infinity)')
ans =
1
>> maple('limit(simplify((a(n))^(1/n)), n=infinity)')
ans =
1
>> maple('limit(n*(1-a(n+1)/a(n)),n=infinity)')
ans =
3
The root and ratio tests tells us nothing, but since the limit is greater than 1, the Raabe criterion tells us that
the series converges.
54
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-9
This is an alternating series. Let us consider the series of moduli and analyze its character:
>> maple ('a: = n - > 1 /(2*n^2+1)');
We apply to this series of positive terms the criteria of comparison of the second kind, comparing with the
convergent series with general term 1/n 2:
>> maple('limit(a(n)/(1/(n^2)),n=infinity)')
ans =
½
As the limit is greater than zero, the given series of positive terms is convergent, so the initial series is absolutely
convergent and, therefore, convergent.
EXERCISE 3-10
55
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
n
is monotone decreasing with limit 0.
Using the Dirichlet test we conclude that the series is convergent.
EXERCISE 3-11
47 49 ×
16 < x < 16
We already know that for values of x in the previous interval the series is convergent. Now we need to analyze the
behavior of the series at the end points of the interval. We first consider x = 49/16:
>> maple ('a: = x - > (4 ^(2*n)) * ((x-3) ^ n) / (n + 2)');
>> maple('simplify(a(49/16))')
56
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
ans =
1/(n+2)
We first apply convergence tests for non-negative series to see if any of them determine the convergence or
divergence of the series.
The ratio, Raabe and root tests do not solve the problem. Next we apply the criterion of comparison of the second
kind, comparing the series of the problem with the divergent harmonic series with general term 1/n:
>> maple('a:=n->1/(n+2)');
>> maple('limit(simplify((a(n))^(1/n)), n=infinity)')
ans =
1
As the limit is greater than zero, the series is divergent.
We now analyze the behavior of the series at the other end point x = 47/16:
>> maple ('a: = x - > (4 ^(2*n)) * ((x-3) ^ n) / (n + 2)');
>> maple ('simplify (a(47/16))')
ans =
(- 1) ^ n / (n + 2)
¥
(-1)n
We have to analyze the alternating series å n + 2 .
n =1
The series with general term (-1)n has bounded partial sums, and the sequence with general term 1 / (n + 2)
is decreasing toward 0. Then, by the Dirichlet test, the alternating series converges. Therefore the interval of
convergence of the power series is the half-closed interval [47/16, 49/16).
EXERCISE 3-12
57
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
58
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
59
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
2 3 4 19 5 49 6 7
2 + 3/2 x + 3/4 x + 7/24 x + 5/48 x + --- x + ---- x + O(x )
480 2880
>> maple('s := multiply(t, v): ')
>> pretty(simple(sym(maple('tpsform(s, x) '))))
2 19 3 4 109 5 6
1 + 3/2 x + 5/4 x + -- x + 7/16 x + --- x + O(x )
24 480
>> maple('s := powexp(x): ')
>> pretty(simple(sym(maple('tpsform(s, x, 7) '))))
2 3 4 5 6 7
1 + x + 1/2 x + 1/6 x + 1/24 x + 1/120 x + 1/720 x + O(x )
>> maple('t := powexp( exp(x) ): ')
>> pretty(simple(sym(maple('tpsform(t, x, 4) '))))
2 3 4
exp (1) + exp (1) x + exp (1) x + 5/6 exp (1) x + or (x)
>> maple('u := powexp( powdiff( powlog(1+x) ) ): ')
>> pretty(simple(sym(maple('tpsform(u, x, 5) '))))
2 3 73 4 5
exp(1)- exp(1) x + 3/2 exp(1) x - 13/6 exp(1) x + -- exp(1) x + O(x)
24
>> maple('powcreate(t(n)=t(n-1)/n,t(0)=0,t(1)=1): ')
>> maple('powcreate(v(n)=v(n-1)/2,v(0)=0,v(1)=1): ')
>> maple('s := reversion(t,v): ')
>> pretty(simple(sym(maple('tpsform(s,x,11) '))))
3 5 7 9 11
x + 1/12 x + 1/80 x + 1/448 x + 1/2304 x + O(x )
>> maple('ts := compose(t,s): ')
>> pretty(simple(sym(maple('tpsform(ts,x) '))))
2 3 4 5 6
x + 1/2 x + 1/4 x + 1/8 x + 1/16 x + O(x )
>> pretty(simple(sym(maple('tpsform(v,x) '))))
2 3 4 5 6
x + 1/2 x + 1/4 x + 1/8 x + 1/16 x + O(x )
>> maple('powcreate(f(n)=f(n-1)/n,f(0)=1): ')
>> maple('powcreate(g(n)=g(n-1)/2,g(0)=0,g(1)=1): ')
>> maple('powcreate(h(n)=h(n-1)/5,h(0)=1): ')
>> maple('k:=evalpow(f^3+g-quotient(h,f)): ')
>> pretty(simple(sym(maple('tpsform(k,x,5) '))))
60
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
61
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
62
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
63
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
64
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-13
Find the Taylor series up to 13th order of sinh (x) at the point x = 0. Also find the Taylor expansion of 1 /(1+x) up to
10th order at the point x = 0. Find the corresponding Laurent and Euler-McLaurin series. Find the Chebychev-Padé
approximation of degree (4,3) in the interval [0,1].
>> pretty(simple(sym(maple('taylor(sinh(x),x=0,13) '))))
65
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
3 5 7 9 11 13
x + 1/6 x + 1/120 x + 1/5040 x + 1/362880 x + 1/39916800 x + O(x )
>> pretty(simple(sym(maple('taylor(1/(1+x),x=0,10) '))))
2 3 4 5 6 7 8 9 10
1 - x + x - x + x - x + x - x + x - x + O(x )
>> maple('with(numapprox): ')
>> pretty(simple(sym(maple('laurent(sinh(x),x=0,13) '))))
3 5 7 9 11 13
x + 1/6 x + 1/120 x + 1/5040 x + 1/362880 x + 1/39916800 x + O(x )
>> pretty(simple(sym(maple('laurent(1/(1+x),x=0,10) '))))
2 3 4 5 6 7 8 9 10
1 - x + x - x + x - x + x - x + x - x + O(x )
Note that the Laurent and Taylor series agree in this case.
>> maple('readlib(eulermac) : ')
>> pretty(simple(sym(maple('eulermac(sinh(x),x,13) '))))
314416268077
------------ cosh(x) - 1/2 sinh(x) + O(cosh(x))
290594304000
>> pretty(simple(sym(maple('eulermac(1/(1+x),x,10) '))))
1 1 1
ln(1 + x) - 1/2 ----- - 1/12 -------- + 1/120 --------
1 + x 2 4
(1 + x) (1 + x)
1 1 1 1
- 1/252 -------- + 1/240 -------- - 1/132 --------- + O(---------)
6 8 10 12
(1 + x) (1 + x) (1 + x) (1 + x)
>> maple(' with(numapprox): ')
>> pretty(simple(sym(maple(' with(orthopoly,T):chebpade(sinh(x),x=0..1,[4,3]) '))))
2
(.03672350156 + .9734745062 x - .02546155828 (2. x - 1.)
3 4
+ .009531362262 (2. x - 1.) - .001730582994 (2. x - 1.) ) /
2
(1.157361890 - .3056450150 x - .009078765811 (2. x - 1.)
3
+ .001929812884 (2. x - 1.) )
66
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
EXERCISE 3-14
Find the Laurent expansion of 1 / (x sin (x)) at the point x = 0 up to degree 10. Also find the Padé approximation of
degree (3,4) at the same point and the Taylor expansion of log (x) at the point x = 2 up to degree 5.
>> pretty(simple(sym(maple('with(numapprox):laurent(1/(x*sin(x)),x=0,10) '))))
-2 2 31 4 127 6 7
x + 1/6 + 7/360 x + ----- x + ------ x + O(x )
15120 604800
>> pretty(simple(sym(maple('with(numapprox):pade(1/(x*sin(x)),x=0,[3,4]) '))))
2
-60 - 3 x
------------
4 2
7 x - 60 x
>> pretty(simple(sym(maple('taylor(log(x),x=2,4) '))))
2 3 4
ln(2) + 1/2 (x - 2) - 1/8 (x - 2) + 1/24 (x - 2) + O((x - 2) )
EXERCISE 3-15
Find the Taylor series at the origin, both in its normal form and in its Poisson form, up to order 5, of the function:
2
f(x, y) = e x + y
67
www.pdfgrip.com
Chapter 3 ■ Numerical Series and Power Series
2 2 2 3 4 2 2 4
1 + x + y + 1/2 x + x y + 1/6 x + 1/2 y + 1/2 y x + 1/24 x
EXERCISE 3-16
Find the Taylor expansion of the function 1 /(2−x) at the point x = 1 up to order 10. Also find the Taylor expansion
of sin(x) at pi/2 up to order 8, and the same for the function log (x) at the point x = 2 up to degree 5.
>> pretty(simple(sym(maple('taylor(1/(2-x),x=1,10)'))))
1 + x - 1 + (x - 1) + (x - 1) + (x - 1) + (x - 1) + (x - 1) + (x - 1) +
8 9 10
(x - 1) + (x - 1) + O((x - 1) )
>> pretty(simple(sym(maple('taylor(sin(x),x=pi/2,8)'))))
2 4 6 8
1 - 1/2(x - 1/2 pi) + 1/24(x - 1/2 pi) - 1/720(x - 1/2 pi) + O((x - 1/2 pi) )
>> pretty(simple(sym(maple('taylor(log(x),x=2,5)'))))
2 3 4 5
log(2) + 1/2 (x - 2) - 1/8 (x - 2) + 1/24 (x - 2) - 1/64 (x - 2) + O((x - 2) )
68
www.pdfgrip.com
Chapter 4
f (a + h ) - f (a )
lim = f ¢(a )
h®0 h
The value of the limit, if it exists, is denoted by f ' (a), and is called the derivative of the function f at the point a.
If f is differentiable at every point of its domain, it is simply said to be differentiable.
The continuity of a function is a necessary condition for its differentiablity, and all differentiable functions are
continuous.
MATLAB can make use of the Maple library is differentiable, whose commands allow you to analyze the
differentiability of a function, and calculate the points where it is not differentiable. The syntax is as follows:
readlib (isdifferentiable): isdifferentiable (expression, var, n, name)
This determines whether the given expression is differentiable to order n in the variable var.
The expression can contain piecewise-defined functions. The optional variable name stores the
sequence of points where the expression is not differentiable.
Exercise 4-1
69
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Figure 4-1.
70
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
You can also determine whether the function is differentiable using the following:
>> maple('readlib(isdifferentiable):')
>> maple('isdifferentiable(piecewise(x=0,0,x*sin(1/x)),x,0)')
ans =
false
Exercise 4-2
Study the differentiability of the function:
x3
f (x) = if x ¹ 0 and f (x) = 0 if x = 0.
|x|
>> maple('f:=x->x^3/abs(x)');
>> maple('limit((f(0+h)- 0)/h,h=0,left)')
ans =
0
>> maple('limit((f(0+h)- 0)/h,h=0,right)')
ans =
0
Thus we see that the derivative at zero exists and has the value zero. The function is plotted in Figure 4-2, from
which we see that it appears to be differentiable at all points of its domain:
>> fplot('x^3/abs(x)',[-1/10,1/10])
71
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Figure 4-2.
72
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
(D@@n)(f)(num) returns the value of the nth derivative of f at the point num
diff (expr, var) returns the derivative of expr with respect to the specified variable
diff (expr, var$ n) returns the derivative of order n of expr with respect to the specified
variable
diff(expr,var1,var2,...,varn):
Returns the partial derivative of expr with respect to the variables var1 to varn in this
order
Diff(expr,var) returns the inert derivative of expr with respect to the variable var
Diff(expr,var1,..,varn) returns the inert partial derivative of expr with respect to var1 to
varn in this order
op (Diff (expr, var)) shows the internal structure of Diff and allows access to its
operands
op (n, Diff(expr,var)) extracts the nth operand (the first is expr, the second is var,
the 3rd, 4th,... are the variables). You can replace any operand with subsop.
convert (expr, D) converts an expression of the form diff (f (x), x) to the form D (f) (x)
convert (expr, diff) converts an expression of the form D (f) (x) to diff (f (x), x)
implicitdiff (equation, vardep, varind) returns the derivative (dy/dx) of the implicit
function specified by the given equation
implicitdiff (equation, vardep, varind, n) returns the derivative of order n (d yn / dxn) of
the implicit function defined by the given equation
implicitdiff(equation,vardep,varind1,...,varindn) returns the derivative (dny/dx1...dxn)
of the implicit function defined by the given equation
implicitdiff (expression, vardep, varindp) returns the derivative of the implicit function
defined by the equation expression = 0
readlib(isdiffentiable):isdifferentiable(expression,var,n): This determines if the
expression is differentiable to order n in the variable var. The expression can contain
piecewise defined functions.
Exercise 4-3
log (sin ( 2 x ) ) , x
tan( x )
,
4 x2 1
3 x2 + 2 (
, log x + x2 + 1 .)
>> pretty(simple(diff('log(sin(2*x))','x')))
2
--------
tan(2 x)
>> pretty(simple(diff('x^tanx','x')))
73
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
tanx
x tanx
------------
x
>> pretty(simple(diff('(4/3)*sqrt((x^2-1)/(x^2+2))','x')))
x
4 -----------------------
2 1/2 2 3/2
(x - 1) (x + 2)
>> pretty(simple(diff('log(x+(x^2+1)^(1/2))','x')))
1
------------
2 1/2
(x + 1)
Exercise 4-4
1 x 2 1+ x
, e ,
x 1- x
>> f='1/x';
>> [diff(f),diff(f,2),diff(f,3),diff(f,4),diff(f,5)]
ans =
-1/x ^ 2 2/x ^ 3 - 6/x ^ 4 24/x ^ 5 - 120/x ^ 6
We see from the pattern here that the nth derivative is given by (-1n)+1n ! .
n
x
>> f='exp(x/2)';
>> [diff(f),diff(f,2),diff(f,3),diff(f,4),diff(f,5)]
ans =
1/2*exp(1/2*x) 1/4*exp(1/2*x) 1/8*exp(1/2*x) 1/16*exp(1/2*x) 1/32*exp(1/2*x)
x /2
Thus the nth derivative is e n .
2
>> f='(1+x)/(1-x)';
>> [simple(diff(f)),simple(diff(f,2)),simple(diff(f,3)),simple(diff(f,4))]
74
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
ans =
2 /(-1+x) ^ 2-4 /(-1+x) ^ 3 12 /(-1+x) ^ 4-48 /(-1+x) ^ 5
2 (n !)
Thus, the nth derivative is given by .
(1 - x )n+1
If f is a function for which f ' (x0) and f '' (x0) both exist, then, if f ' (x0) = 0 and f '' (x0) < 0, the function f has a local
maximum at the point (x0, f (x0)).
If f is a function for which f ' (x0) and f '' (x0) both exist, then, if f ' (x0) = 0 and f '' (x0) > 0, the function f has a local
minimum at the point (x0, f (x0)).
If f is a function for which f ' (x0), f '' (x0) and f ''' (x0) exist, then, if f ' (x0) = 0 and f '' (x0) = 0 and f ''' (x0) ≠ 0, the
function f has a turning point at the point (x0, f(x0)).
If f is differentiable, then the values of x for which the function f is increasing are those for which f ' (x) is greater
than zero.
If f is differentiable, then the values of x for which the function f is decreasing are those for which f ' (x) is less than zero.
If f is twice differentiable, then the values of x for which the function f is concave are those for which f '' (x) is
greater than zero.
If f is twice differentiable, then the values of x for which the function f is convex are those for which f '' (x) is less
than zero.
Exercise 4-5
ans =
-12
>> subs(f,-1)
ans =
20
Thus the slope of the tangent line at the point x = −1 is − 12, and the function at x = −1 has the value 20.
Therefore the equation of the tangent to the curve at the point (− 1, 20) is:
y - 20 = - 12 ( x - (- 1))
We plot the curve and its tangent on the same axes (see Figure 4-3):
>> fplot('[2*x^3+3*x^2-12*x+7, 20-12*(x - (-1))]',[-4,4])
Figure 4-3.
76
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
To calculate the horizontal tangent to the curve y = f(x), we find the values x0 for which f ' (x0) = 0. The equation of
this tangent will then be y = f(x0).
To calculate the vertical tangent to the curve y = f(x), we find the values x0 for which f ' (x0) = ¥. The equation of
this tangent will then be x = x0.
>> g ='(x^2-x+4) /(x-1)'
>> solve(diff(g))
ans =
[3]
[-1]
>> subs(g,3)
ans =
5
>> subs(g,-1)
ans =
-3
The two horizontal tangents will have equations:
y = g '[−1] (x + 1) −3, that is, y = −3
y = g '[3] (x − 3) + 5, that is, y = 5
The horizontal tangents are not asymptotes because the corresponding values of x0 are finite (they are - 1 and 3).
We now find the vertical tangents. To do this, we calculate the values of x for which g'(x) = ¥ (i.e. values for which
the denominator of g ' is zero, but does not cancel with the numerator):
>> solve('x-1')
ans =
1
Therefore, the vertical tangent has equation x = 1.
For x = 1, the value of g(x) is infinite (see below), so the vertical tangent is a vertical asymptote.
subs(g,1)
Error, division by zero
77
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Figure 4-4.
78
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Exercise 4-6
Find the asymptotes, maxima, minima, inflection points, intervals of growth and decrease and intervals of
concavity and convexity for the function:
x3
f (x) =
x -1
2
>> f='x^3/(x^2-1)'
f =
x^3/(x^2-1)
>> syms x, limit(x^3/(x^2-1),x,inf)
ans =
NaN
Therefore, there are no horizontal asymptotes. To see if there are any vertical asymptotes, let us look at the
values of x that make y infinite:
>> solve('x^2-1')
ans =
[1]
[-1]
Thus the vertical asymptotes are the lines x = 1 and x = −1. Now let us see if there are any oblique asymptotes:
>> limit(x^3/(x^2-1)/x,x,inf)
ans =
1
>> limit(x^3/(x^2-1)-x,x,inf)
ans =
0
The line y = x is an oblique asymptote. Now we shall find the maxima and minima, inflection points and intervals
of concavity:
>> solve (diff (f))
79
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
ans =
[ 0]
[ 0]
[3 ^(1/2)]
[^(1/2) - 3]
The first derivative vanishes at x = 0, x = ÷3 and x = −÷3. These include maximum and minimum candidates.
To verify if they are maxima or minima, we find the value of the second derivative at these points:
>> [numeric(subs(diff(f,2),0)),numeric(subs(diff(f,2),sqrt(3))),
numeric(subs(diff(f,2),-sqrt(3)))]
ans =
0 2.5981 - 2.5981
Therefore, at x = –÷3 there is a maximum and at x =÷3 there is a minimum. At x = 0 we know nothing:
>> [numeric (subs (f, sqrt (3))), numeric (subs (f, - sqrt (3)))]
ans =
2.5981 - 2.5981
Therefore, the maximum point is (−÷3, −2.5981) and the minimum point is (Ö3, 2.5981).
We will now analyze the points of inflection:
>> solve(diff(f,2))
[ 0]
[ i*3^(1/2)]
[-i*3^(1/2)]
The only possible turning point occurs at x = 0, and as f (0) = 0, the possible turning point is (0,0):
>> subs (diff(f,3), 0)
ans =
-6
As the third derivative at x = 0 is non-zero, the origin really is a turning point:
>> pretty(simple(diff(f)))
2 2
x (x - 3)
------------
2 2
(x - 1)
80
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
The curve is increasing when y' > 0, i.e., in the intervals (-¥,-÷ 3) and (÷ 3,¥).
The curve is decreasing when y ' < 0, i.e., in the intervals
(-÷ 3,-1), (-1,0), (0,1) and (1, ÷ 3).
>> pretty(simple(diff(f,2)))
2
x (x + 3)
2 ------------
2 3
(x - 1)
The curve will be concave when y" > 0, i.e., in the intervals (-1,0) and (1, ¥ ).
The curve is convex when y" < 0, i.e. in the intervals (0,1) and (- ¥ , - 1).
The curve has a horizontal tangent at the three points at which the first derivative is zero. The equations of the
horizontal tangents are y = 0, y = -2.5981 and y = 2.5981.
The curve has vertical tangents at the points that make the first deriviative infinite. These are x = 1 and x = −1.
Therefore, the vertical tangents coincide with the two vertical asymptotes.
We plot the curve together with its asymptotes (see Figure 4-5):
>> fplot('[x^3/(x^2-1),x]',[-5,5,-5,5])
Figure 4-5.
81
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
We can also represent the curve, its asymptotes and its horizontal and vertical tangents in the same graph
(Figure 4-6):
>> fplot('[x^3/(x^2-1),x,2.5981,-2.5981]',[-5,5,-5,5])
Figure 4-6.
Exercise 4-7
Write a positive number a as the sum of two summands in such a way that the sum of their cubes is minimal.
Let x be one of the summands. The other will be a − x. We require the sum x 3 + (a − x)3 to be minimized.
>> f='x^3+(a-x)^3'
f =
x^3+(a-x)^3
82
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
>> solve(diff(f,'x'))
ans =
1/2 * a
The possible maximum or minimum is at x = a/2. We use the second derivative to verify that it is indeed a minimum:
>> subs(diff(f,'x',2),'a/2')
ans =
3 * a
As a > 0 (by hypothesis), 3a > 0, which means x = a/2 is a minimum.
Therefore the two summands must be equal to a/2.
Exercise 4-8
Suppose you want to purchase a rectangular plot of 1600 square meters and then fence it. Knowing that the
fence costs 200 cents per meter, what dimensions must the plot of land have to ensure that the fencing is most
economical?
If the surface area is 1600 square feet and one of its dimensions, unknown, is x, and the other will be 1600/x.
The perimeter of the rectangle is p (x) = 2 x + 2 (1600/x), and the cost is given by f(x) = 200 p(x):
>> f ='200 * (2 * x + 2 *(1600/x))'
f =
200 * (2 * x + 2 *(1600/x))
This is the function to minimize:
>> solve (diff (f))
ans =
[40]
[-40]
The possible maxima and minima are presented for x = − 40 and x = 40. To determine their nature, we look at
the second derivative at these points.
>> [subs (diff(f,2), 40), subs (diff(f,2), - 40)]
ans =
20 -20
83
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
x = 40 is a minimum, and x = − 40 is a maximum. Thus, one of the sides of the rectangular field must be 40
meters, and the other will measure 1600/40 = 40 meters. Therefore the optimal rectangle is a square with sides
of 40 meters.
Exercise 4-9
84
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
>> pretty(sym(maple('D[2,2](f)')))
2 2 2 2 2
(x,y) -> - sin(x y) x - 4 cos(x y ) x y - 2 sin(x y ) x
>> pretty(sym(maple('D[1,2](f)')))
2 3 2
(x,y) -> - sin(x y) y x + cos(x y) - 2 cos(x y ) y x - 2 sin(x y ) y
>> pretty(sym(maple('D[2,1](f)')))
2 3 2
(x,y) -> - sin(x y) y x + cos(x y) - 2 cos(x y ) y x - 2 sin(x y ) y
>> pretty(sym(maple('D[1,1,2,2](f)')))
2 2 2 6 2
(x,y) -> sin(x y) y x - 4 cos(x y) y x - 2 sin(x y) + 4 cos(x y ) y x
2 4 2 2
+ 18 sin(x y ) y x - 12 cos(x y ) y
Exercise 4-10
Calculate the implicit derivative of y with respect to x, where x and y satisfy the following equation:
x3 + 3 x2y3 + 5xy + 6y3 = 8
First, we will use the operator D to differentiate the equation:
>> pretty (sym (maple ('D(x^3+3*x^2*y^3+5*x*y+6*y^3=8)')))
2 3 2 2 2
3 D(x) x + 6 D(x) x y + 9 x D(y) y + 5 D(x) y + 5 x D(y) + 18 D(y) y = 0
Now we replace D (x) with 1 and solve for D (y):
>> pretty (sym (maple ('subs (D (x) = 1, D(x^3+3*x^2*y^3+5*x*y+6*y^3=8))')))
2 3 2 2 2
3 x + 6 x y + 9 x D(y) y + 5 y + 5 x D(y) + 18 D(y) y = 0
85
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Exercise 4-11
86
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
87
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Figure 4-7.
MATLAB has the command maple(‘completesquare’), implemented in the student library, which allows you to
complete the square of an expression with respect to a specified variable.
We complete the square of the defining expression of the ellipse with respect to the variable y as follows:
>> pretty(maple('completesquare(2*x^2-2*x*y+y^2+x+2*y+1,y)'))
2 2
(y - x + 1) + x + 3 x
Now, it is trivial to separate the variable y for the given equation: y = x - 1 + - x 2 - 3 x . But we have an
expression within a square root, which we rewrite as follows:
>> pretty(maple('completesquare(-x^2-3*x,x)'))
2
- (x + 3/2) + 9/4
We obtain the equation of the curve
2
æ 3ö
y = x -1+ -ç x + ÷ + 9 / 4
è 2ø
3
which is easily parametrized by putting x + 3/2 = 3/2 sin (t ), so that x = * ( sin (t ) - 1) , y = 3 / 2 * (sin (t ) + cos (t )) - 5 / 2.
2
Now we can create the graph of the curve, the asymptotes and the coordinate axes (the OY axis is represented in
parametric form as x = 0, y = t ).
88
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
>> fplot('[x+3/2-1,x+3/2-4,0]',[-4,2,-6,2]);
>> hold on;
>> t=(0:.1:2*pi);
>> x=3/2*(sin(t)-1);
>> y=3/2*(sin(t)+cos(t))-5/2;
>> plot(x,y);
>> t=(-6:.1:2);
>> x=zeros(size(t));
>> y=t;
>> plot(x,y)
¶f f (a + h , b ) - f (a , b )
(a , b) = lim
¶x h®0 h
Similarly, the partial derivative of f with respect to the variable y at the point (a, b) is defined as:
¶f f (a , b + h ) - f (a , b )
(a , b) = lim
¶y h®0 h
Generally speaking, we can define the partial derivative with respect to any variable for a function of n variables.
Given the function f: Rn→ R, the partial derivative of f with respect to the variable xi (i = 1, 2,..., n) at the point
(a1, a2, ..., an) is defined as follows:
The function f is differentiable if all partial derivatives with respect to xi (i = 1, 2,..., n) exist and are continuous.
All differentiable functions are continuous, and if a function is not continuous, it cannot be differentiable.
The directional derivative of the function f with respect to the vector v = (v1,v2,...,vn) is defined as the following
scalar product:
æ¶f ¶f ¶f ö
(Df )v = ç , ,..., ÷ · (v1 ,v2 ,...,vn ) = ( Ñf ) · v
¶
è 1 x ¶ x 2 ¶ xn ø
æ¶f ¶f ¶f ö
Ñf = ç , ,..., ÷ is called the gradient vector of f.
è ¶ x1 ¶ x 2 ¶ xn ø
The directional derivative of the function f with respect to the vector v = (dx1,dx2,...,dxn) is called the total
differential of f. Its value is:
æ¶f ¶f ¶f ö
Df = ç dx1 + dx 2 + ... + dx n ÷
¶
è 1 x ¶ x 2 ¶ x n ø
89
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
At the beginning of the chapter (under the heading “Calculating derivatives”) we defined the functions that
MATLAB implements to find partial derivatives. Using these functions, we have the following:
∂f
diff (f (x, y, z...), x) = —— = maple ('D [1] (f)(x,y,z,...)') = maple ('diff (f(x,y,z,...), x)')
∂x
∂nf n)
diff (f(x,y,...), x, n) = —— = maple('D[1,1,...1] (f)(x,y,...)') = maple ('diff (f(x,y,...), x$n)')
∂xn
∂f
diff(f(x1,x2,...),xj) = ———— maple('D[ j ](f)(x1,x2,...)') = maple(‘diff(f(x1,x2,..), xj)’)
∂xj
∂nf n)
diff(f(x1,...),xj,n)= ———— = maple('D[j, j,...j](f)(x1,...)') = maple(‘diff(f(x1,...), xj $n)’)
∂xjn
∂∂∂
maple ('diff (f(x,y,z,...), x, y, z,...)') = ———————— ...f = maple('D[1,2,3..])(f)(x,y,z,...)')
∂x ∂y ∂z
∂2f
maple ('diff (f(x,y), x, y)') = ——————— = maple ('D [1,2] (f)(x,y)')
∂x∂y
∂(k)f
maple(‘diff(f(x1,...),xn,xm,xp)’)= ————————— = maple(‘D[n,m,p,...](f)(x1,x2,...)’)
∂xn∂xm∂xp ... (k = number of variables)
Exercise 4-12
90
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
ans =
0
The limits of the two previous expressions as h 0 and k 0, respectively, are both zero.
We see that the two partial derivatives at the origin are the same and have value zero. But the function is not
differentiable at the origin, because it is not continuous at (0,0), since there is no limit as (x, y ) tends to (0,0 ):
>> maple('limit((m*x)^2/(x^2+(m*x)^2),x=0)')
ans =
m ^ 2 /(1+m^2)
The limit does not exist at (0,0), because approaching the origin along different straight lines y = mx yields
different results depending on the parameter m.
Exercise 4-13
The function is differentiable if it has continuous partial derivatives at every point. We will consider any point other
than the origin and calculate the partial derivative with respect to the variable x :
>> maple('f:=(x,y)->(2*x*y)/(x^2+y^2)^(1/2)');
>> pretty(simple(sym(maple('D[1](f)(x,y)'))))
3
y
2 --------------
2 2 3/2
(x + y )
Now, let’s see if this partial derivative is continuous at the origin:
>> maple('limit(2*(m*x)^3/(x^2+(m*x)^2)^(3/2),x=0)')
ans =
2*m^3/(1+m^2)^(3/2)
The limit does not exist at (0,0), because approaching the origin along different straight lines y = mx yields
different results depending on the parameter m. Therefore, the partial derivative is not continuous at the origin.
We conclude that the function is not differentiable.
91
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
However, the function is continuous, since the only problematic point is the origin, and the limit of the function at
the origin is 0 = f (0,0):
>> maple('limit(limit(f(x,y),x=0),y=0)')
ans =
0
>> maple('limit(limit(f(x,y),y=0),x=0)')
ans =
0
>> maple('limit(f(x,(m*x)),x=0)')
ans =
0
>> maple('limit(f(x,(m*x)^(1/2)),x=0)')
ans =
0
>> maple('limit(f((r*cos(a)),(r*sin(a))),r=0)')
ans =
0
The iterated limits and the directional limits are all zero, and by changing the function to polar coordinates, the
limit at the origin turns out to be zero, which coincides with the value of the function at the origin.
This is therefore an example of a non-differentiable continuous function.
Exercise 4-14
92
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Grouping the terms together, we obtain the expression of the total differential:
>> pretty(sym(maple('collect(D(f(x,y)),D(x))')))
2 2
(3 x - 4 x y) D(x) - 2 x D(y)
>> maple('g:=(x,y)->exp(a*x)*cos(b*y)');
>> pretty(sym(maple('dotprod([D[1](g)(x,y),D[2](g)(x,y)],[dx,dy])')))
a exp(a x) cos(b y) dx - exp(a x) sin(b y) b dy
Exercise 4-15
1
then ¶ f2 + ¶ f2 + ¶ f2 = 0.
2 2 2
Verify that if f ( x , y , z ) =
x +y +z
2 2 2
¶x ¶y ¶z
>> f='1/(x^2+y^2+z^2)^(1/2)';
>> simplify(symop(diff(f,'x',2),'+',diff(f,'y',2),'+',diff(f,'z',2)))
ans =
0
Another path would be the following:
>> maple('f:=x->1/(x^2+y^2+z^2)^(1/2)');
>> simplify(sym('diff(f(x,y,z),x$2)+diff(f(x,y,z),y$2)+diff(f(x,y,z),z$2)'))
ans =
0
Alternatively, another way to do the same is as follows:
>> maple('simplify(D[1,1](f)(x,y,z)+D[2,2](f)(x,y,z)+D[3,3](f)(x,y,z))')
ans =
0
Exercise 4-16
1
f (x , y ,z) =
x + y2 + z2
2
>> maple('f:=(x,y,z)->1/(x^2+y^2+z^2)^(1/2)');
>> maple('dotprod([D[1](f)(2,1,1),D[2](f)(2,1,1),D[3](f)(2,1,1)],[1,1,0])')
ans =
-1/12*6^(1/2)
Exercise 4-17
calculate:
¶ f ¶ f ¶2f ¶2f ¶2f ¶3f ¶4f ¶5f
, , , , , , and
¶x ¶y ¶x¶y ¶x ¶y ¶x¶y ¶x ¶y
2 2 2 2 2
¶ x3 ¶ y2
ans =
0.6745
>> pretty(factor(sym('D[1,2](f)(x,y)')))
2 2
1/16 exp(- 1/8 x - 1/8 y )
2 2
(y x cos(x) + y x sin(y) + 8 y cos(x) sin(x) - 8 x sin(y) cos(y))
>> pretty(factor(sym('D[1,2](f)(pi/3,pi/6)')))
2 1/2
1/576 pi exp(- 5/288 pi ) (pi - 12 3 )
>> numeric('D[1,2](f)(pi/3,pi/6)')
ans =
-0.0811
>> pretty(factor(sym('D[2,1](f)(x,y)')))
2 2
1/16 exp(- 1/8 x - 1/8 y )
2 2
(y x cos(x) + y x sin(y) + 8 y cos(x) sin(x) - 8 x sin(y) cos(y))
>> pretty(factor(sym('D[2,1](f)(pi/3,pi/6)')))
2 1/2
1/576 pi exp(- 5/288 pi ) (pi - 12 3 )
>> numeric('D[2,1](f)(pi/3,pi/6)')
ans =
-0.0811
>> pretty(factor(sym('D[1,1](f)(x,y)')))
2 2 2 2 2 2 2 2
1/16 exp(- 1/8 x - 1/8 y ) (- 36 cos(x) - 4 sin(y) + x cos(x) + x sin(y)
2
+ 16 x cos(x) sin(x) + 32 sin(x) )
95
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
>> pretty(factor(sym('D[1,1](f)(pi/3,pi/6)')))
2 2 1/2
1/288 exp(- 5/288 pi ) (252 + pi + 24 pi 3 )
>> numeric('D[1,1](f)(pi/3,pi/6)')
ans =
1.1481
>> pretty(factor(sym('D[2,2](f)(x,y)')))
2 2 2 2 2 2 2 2
1/16 exp(- 1/8 x - 1/8 y ) (- 4 cos(x) - 36 sin(y) + y cos(x) + y sin(y)
2
- 16 y sin(y) cos(y) + 32 cos(y) )
>> pretty(factor(sym('D[2,2](f)(pi/3,pi/6)')))
2 2 1/2
- 1/1152 exp(- 5/288 pi ) (- 1008 - pi + 48 pi 3 )
>> numeric('D[2,2](f)(pi/3,pi/6)')
ans =
0.5534
>> pretty(factor(sym('D[1,1,1,2](f)(x,y)')))
2 2 2 2
1/256 exp(- 1/8 x - 1/8 y ) (- 108 y x cos(x) - 12 y x sin(y)
3 2 3 2
- 608 y cos(x) sin(x) + y x cos(x) + y x sin(y)
2 2
+ 24 y x cos(x) sin(x) + 96 y x sin(x) + 96 x sin(y) cos(y)
3
- 8 x sin(y) cos(y))
>> pretty(factor(sym('D[1,1,1,2](f)(pi/3,pi/6)')))
2 3 1/2 2 1/2
1/82944 pi exp(- 5/288 pi ) (pi + 12 3 pi + 756 pi - 5616 3 )
>> numeric('D[1,1,1,2](f)(pi/3,pi/6)')
96
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
ans =
-0.2271
>> pretty(factor(sym('D[1,1,2,2](f)(x,y)')))
2 2 2 2 2 2
1/256 exp(- 1/8 x - 1/8 y ) (144 cos(x) + 144 sin(y) - 4 x cos(x)
2 2 2 2 2
- 36 x sin(y) - 64 x cos(x) sin(x) - 128 sin(x) - 36 y cos(x)
2 2 2 2 2 2 2 2 2
- 4 y sin(y) + y x cos(x) + y x sin(y) + 16 y x cos(x) sin(x)
2 2 2 2
+ 32 y sin(x) + 64 y sin(y) cos(y) - 16 y x sin(y) cos(y) - 128 cos(y)
2 2
+ 32 x cos(y) )
>> pretty(factor(sym('D[1,1,2,2](f)(pi/3,pi/6)')))
1/165888*
2 2 1/2 4 3 1/2
exp(- 5/288 pi ) (- 77760 + 1260 pi - 1728 pi 3 + pi - 24 pi 3 )
>> numeric('D[1,1,2,2](f)(pi/3,pi/6)')
ans =
-0.3856
>> pretty(factor(sym('D[1,1,1,2,2](f)(x,y)')))
2 2 2
- 1/1024 exp(- 1/8 x - 1/8 y ) (2432 cos(x) sin(x) + 432 x cos(x)
2 2 3 2 3 2
+ 432 x sin(y) - 384 x sin(x) - 4 x cos(x) - 36 x sin(y)
2 2 2 2
+ 192 y x sin(y) cos(y) + 96 y x sin(x) - 108 y x cos(x)
2 3 2 2 3 2 2 2 2 2
+ y x cos(x) + y x sin(y) + 24 y x cos(x) sin(x) - 12 y x sin(y)
2 2 3 2
- 608 y cos(x) sin(x) - 96 x cos(x) sin(x) + 32 x cos(y)
2 3
- 384 x cos(y) - 16 y x sin(y) cos(y))
97
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
>> pretty(factor(sym('D[1,1,1,2,2](f)(pi/3,pi/6)')))
2
1/1990656 exp(- 5/288 pi )*
1/2 3 2 1/2 4 1/2 5
(- 1181952 3 + 233280 pi - 1764 pi + 8208 pi 3 + 12 pi 3 - pi )
>> numeric('D[1,1,1,2,2](f)(pi/3,pi/6)')
ans =
-0.5193
The advantage of using the operator D instead of the command diff is that it is much faster at directly calculating
the value of the derivative at a point.
First, suppose that the determinant of H is non-zero at the point (a1, a2, ..., an). In this case, we say that the point is
non-degenerate and, in addition, we can determine the nature of the extreme point via the following conditions:
If the Hessian matrix at the point (a1, a2, ..., an) is positive definite, then the function has a minimum at that point.
If the Hessian matrix at the point (a1, a2, ..., an) is negative definite, then the function has a maximum at that point.
In any other case, the function has a saddle point at (a1, a2, ..., an).
If the determinant of H is zero at the point (a1, a2, ..., an), we say that the point is degenerate.
MATLAB includes a specific function which calculates the determinant of the Hessian of a function of several
variables. Its syntax is:
maple ('hessian(function, [x1,x2,..., xn])')
98
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
MATLAB also includes a function which calculates the gradient of a function of several variables. Its syntax is as
follows:
maple ('grad(function, [x1,x2,..., xn])')
Exercise 4-18
f ( x , y ) = -120 x 3 - 30 x 4 + 18 x 5 + 5 x 6 + 30 xy 2 .
First, we find the extreme points, setting the partial derivatives (the components of the gradient vector of f ) to
zero and solving the resulting system:
>> maple('f:=(x,y)->-120*x^3-30*x^4+18*x^5+5*x^6+30*x* y^2 ');
>> maple('solve({diff(f(x,y),x)=0,diff(f(x,y),y)=0},{x,y})')
ans =
{x = 0, y = 0}, {y = 0, x = 2}, {y = 0, x = -3}, {y = 0, x = -2}
Thus, the extreme points are: (− 2,0), (2,0), (0,0) and (− 3.0).
Alternatively, we can find the gradient vector and solve the system obtained by setting it equal to zero:
>> pretty(sym(maple('grad(f(x,y),[x,y])')))
[-360*x^2-120*x^3+90*x^4+30*x^5+30*y^2, 60*x*y]
>> maple('solve({-360*x^2-120*x^3+90*x^4+30*x^5+30*y^2=0, 60*x*y=0},{x,y})')
ans =
{x = 0, y = 0}, {y = 0, x = 2}, {y = 0, x = -3}, {y = 0, x = -2}
To classify the extreme points, we construct the Hessian matrix, and calculate its value at each point:
>> pretty(sym(maple('hessian(-120*x^3-30*x^4+18*x^5+5*x^6+30*x* y^2,[x,y])')))
[-720*x-360*x^2+360*x^3+150*x^4, 60*y]
[ 60*y, 60*x]
You can also find the Hessian matrix in the following way:
>> pretty(sym(maple('array([[D[1,1](f)(x,y),D[1,2](f)(x,y)],[D[1,2](f)(x,y),
D[2,2](f)(x,y)]])')))
[-720*x-360*x^2+360*x^3+150*x^4, 60*y]
[ 60*y, 60*x]
99
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
But what we need is the Hessian matrix M as a function of x and y so we can evaluate it at certain points:
>> maple ('m:=(x,y) - >
hessian(-120*x^3-30*x^4+18*x^5+5*x^6+30*x* y^2,[x,y])');
>> pretty (sym (maple ('subs (x = 0, y = 0, M(x,y))')))
[0, 0]
[0, 0]
The origin turns out to be a degenerate point, as the determinant of the Hessian matrix is zero at (0,0).
>> pretty (sym (maple ('subs (x =-2, y = 0, M(x,y))')))
[- 480, 0]
[0, - 120]
>> maple ('H: = subs (x =-2, y = 0, M(x,y))');
>> maple('det(H)')
ans =
57600
>> maple('eigenvals(H)')
ans =
-480, - 120
The Hessian matrix at the point (− 2,0) has non-zero determinant, and is also negative definite, because all its
eigenvalues are negative. Therefore, the point (− 2,0) is a maximum of the function.
>> H = sym (maple ('subs (x = 2, y = 0, M(x,y))'))
H =
[2400, 0]
[ 0, 120]
>> determ(H)
ans =
288000
>> eigensys(H)
ans =
[120].
[2400].
100
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
The Hessian matrix at the point (2,0) has non-zero determinant, and is furthermore positive definite, because all
its eigenvalues are positive. Therefore, the point (2,0) is a minimum of the function.
>> H = sym (maple ('subs (x = 3, y = 0, M(x,y))'))
H =
[1350, 0]
[ 0, -180]
>> determ(H)
ans =
-243000
>> eigensys(H)
ans =
[-180.]
[1350]
The Hessian matrix at the point (− 3,0) has non-zero determinant, and, in addition, is neither positive definite nor
negative definite, because its eigenvalues are not all positive or all negative. Therefore (− 3.0) is a saddle point of
the function.
Exercise 4-19
f ( x , y , z ) = x 2 + xy + y 2 + z 2
First, we set the partial derivatives (the components of the gradient vector of f) to zero and solve the resulting
system:
>> maple ('f:=(x,y) - > x ^ 2 + y ^ 2 + z ^ 2 + x * y ');
>> maple ('solve ({diff (f (x, y, z), x) = 0, diff (f (x, y, z), y) = 0, diff (f (x, y, z),
z)}, {x, y, z})')
ans =
{z = 0, y = 0, x = 0}
The single extreme point is the origin (0,0,0). We determine what kind of extreme point it is. To do this, we
calculate the Hessian matrix and express it as a function of x, y and z :
>> maple ('h:=(x,y,z) - > hessian (f (x, y, z), [x, y, z])');
>> M = maple ('H (x, y, z)')
101
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
M =
[2, 1, 0]
[1, 2, 0]
[0, 0, 2]
>> determ(M)
ans =
6
We see that the Hessian matrix is constant (it does not depend on the point (x,y,z )), therefore its value at the
origin is already found. The determinant is non-zero, so the origin is non-degenerate.
>> eigensys (M)
ans =
[1]
[2]
[3]
The Hessian matrix at the origin is positive definite, because all its eigenvalues are positive. Thus we conclude
that the origin is a minimum of the function.
The extreme points are found by solving the system by setting the components of the gradient vector of L to zero,
that is, Ñ L(x1,x2,...,xn,l) = (0,0,...,0). Which translates into:
æ¶L ¶L ¶L ¶L ¶L ¶Lö
ÑL =ç , ,..., , , ,..., ÷ = ( 0 , 0 ,......,0 )
¶
è 1 x ¶ x 2 ¶ x n ¶ l1 ¶ l2 ¶ ln ø
102
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
By setting the partial derivatives to zero and solving the resulting system, we obtain the values of x1, x2,..., xn, l1,
l2,...,lk corresponding to possible maxima and minima.
To determine the nature of the points (x1, x2,..., xn) found above, the following bordered Hessian matrix is used:
The nature of extreme points can be determined by studying the set of bordered Hessian matrices:
é ¶2 f ¶2f ¶ gi ù
ê
é¶ f ¶ gi ù ê ¶ x1
2
¶ x1 ¶ x 2 ¶ x1 úú
ê¶ x2 ¶ x1 ú ê ¶2f ¶2 f ¶ g i ú ...
H1 = ê 1 ú H2 = ê ú Hn = H
ê ¶ gi ú ê ¶ x1 ¶ x 2 ¶ x 22 ¶ x2 ú
ê¶ x 0 ú
ë 1 û ê ¶g ¶ gi ú
ê i
0 ú
êë ¶ x1 ¶ x2 úû
For a single restriction g1, if H1 < 0, H2 < 0, H3 < 0,..., H < 0, then the extreme point is a minimum.
For a single restriction g1, if H1 > 0, H2 < 0, H3 > 0, H4 < 0, H5 > 0, ... then the extreme point is a maximum.
For a collection of restrictions gi(x1,..., xn) (i = 1, 2,..., k) the lower right 0 will be a block of zeros and the conditions
for a minimum will all have sign (-1)k, while the conditions for a maximum will have alternating signs with H1 having
sign (-1)k+1. When considering several restrictions at the same time, it is easier to determine the nature of the extreme
point by simple inspection.
MATLAB also offers some commands for function optimization, all of which require the prior use of the
command maple. Among them are the following:
maximize (expression) is the maximum value of the given expression with respect to
all of its variables
maximize(expression,{var1,...,varn}) is the maximum value of the given expression
with respect to the specified variables
minimize (expression) is the minimum value of the given expression with respect to all
of its variables
minimize(expression,{var1,...,varn}) is the minimum value of the given expression with
respect to the specified variables
minimize(expression,{var1,...,varn},{a1..b1,..., an...bn}) is the minimum value of the
given expression with respect to the specified variables, limiting the ranges of the
variables var1,..., varn to the intervals [a1, b1],..., [an, bn] (by default the variable
ranges are infinite)
103
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
In addition, MATLAB has the command extrema, which gives us the extreme points of a function (the objective
function) subject to a few restrictions, determined by the method of Lagrange multipliers. Its syntax is as follows:
maple ('extrema(fobjective, {restrictions}, {variables},s)')
The variable s contains the list of candidate extreme points. To use this function one must read it from a Maple
library using the command maple (‘readlib(extrema)’). This library also includes the following commands (all
require the prior use of the command maple):
readlib (extrema) loads into the memory a library of commands for optimization with
restrictions
extrema(expression,equation,var) returns candidate extreme points for the function
defined by the given expression in the specified variable under the constraint defined
by the equation, using the method of Lagrange multipliers
extrema(expression,equation) finds candidate extreme points of the expression for all
variables subject to the restriction defined by the equation
extrema(expression, variable, {}) finds candidate extreme points for the expression in
the given variable and without restrictions
extrema(expression,{equation1,...,equationn},{variable1,...,variablen}) finds candidate
extreme points for the given expression with respect to the specified variables subject
to the constraints given by equation1,..., equationn
extrema(expression,{equ1,...,equn},{variable1,...,variablen},n) The variable n returns
the set of values obtained by evaluating the expression (objective function) at every
candidate extreme point. These values will be used to differentiate between minima
and maxima.
Exercise 4-20
ans =
-2 ^(1/2)
Then, at the point (−Ö2 /2, 0,−Ö2 /2) the function has a minimum, and there is a maximum at the point
(Ö2/2,0,Ö2/2).
Exercise 4-21
f (x , y ,z) = x 2 + y 2 - z
105
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
Thus the extreme points are (Ö8, Ö8, −2Ö8 + 10) and (−Ö8, −Ö8, 2Ö8 + 10). Now let’s determine the nature of
these points. To this end, we evaluate the objective function at these points:
>> maple ('obj:=(x,y,z) - >(x^2+y^2) ^(1/2)-z');
>> numeric(maple('obj(-sqrt(8),-sqrt(8),sqrt(8)+10)'))
ans =
-8.8284
>> numeric(maple('obj(sqrt(8),sqrt(8),-sqrt(8)+10)'))
ans =
-3.1716
Thus, the function has a minimum at the point (−Ö8, −Ö8, 2Ö8 + 10) and a maximum at the point (Ö8, Ö8, − 2Ö8 + 10).
Exercise 4-22
Find the dimensions of the rectangular cuboid of maximum volume which has a surface area of 10 square meters.
If x, y, z are the dimensions of the block, its volume will be V = xyz. As we know that it has a surface area of
10 square meters, the restriction will be 2xy + 2xz + 2yz = 10. We therefore have to maximize the objective
function V = xyz subject to the condition 2xy + 2xz + 2yz − 10 = 0. We will use the method of Lagrange multipliers:
>> maple('func:=x*y*z');
>> maple('g1:=2*x*y+2*x*z+2*y*z-10');
>> maple ('readlib (extrema)');
>> maple ('allvalues (extrema(func,{g1},{x,y,z}, a))')
ans =
{5/9 * 15 ^(1/2)}, {- 5/9 * 15 ^(1/2)}
>> maple ('a')
ans =
{{z = 1/3 * 15 ^(1/2), y = 1/3 * 15 ^(1/2), x = 1/3 * 15 ^(1/2)}},
{{z = - 1/3 * 15 ^(1/2), y = - 1/3 * 15 ^ (1/2) x = - 1/3 * 15 ^(1/2)}}
106
www.pdfgrip.com
Chapter 4 ■ Derivatives and Applications. One and Several Variables
107
www.pdfgrip.com
Chapter 5
é ¶ F1 ¶ F1 ¶ F1 ù
ê¶ x .......
¶ x2 ¶ xn ú
ê 1 ú
ê ¶ F2 ¶ F2 ¶ F2 ú
ê ....... ¶ ( F1 , F2 ,..., Fn )
J = ê ¶ x1 ¶ x2 ¶ xn úú =
¶ ( x1 , x 2 ,..., x n )
ê........ ....... ....... ....... ú
ê ú
ê ¶ Fn ¶ Fn
.......
¶ Fn ú
êë ¶ x1 ¶ x2 ¶ xn úû
The Jacobian of a vector function is an extension of the concept of a partial derivative for a single-component
function.
MATLAB has the command jacobian which enables you to calculate the Jacobian matrix of a function:
maple (‘jacobian([F1,F2,...,Fn],[x1,x2,...,xn])’) returns J = ∂ ( F1 , F2 ,..., Fn )
∂ ( x1 , x 2 ,..., x n )
109
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-1
ans =
[exp (x), 0, 0]
[0,-sin (y), 0]
[0, 0, cos (z)]
>> maple ('m:=(x,y,z) - > array ([[exp(x), 0, 0], [0, -sin(y), 0], [0, 0, cos(z)]])');
>> sym(maple('m(0,-Pi/2,0)'))
ans =
[1, 0, 0]
[0, 1, 0]
[0, 0, 1]
We see that the Jacobian matrix evaluated at the given point is the identity matrix.
g : U Ì R n ® R m and f : V Ì R m ® R p .
where U and V are open and consider the composite function f g : R n ® R p .
If g is differentiable at x 0 and f is differentiable at y 0 = g ( x 0 ) , then f g is differentiable at x 0 and we have the
following:
D( f g )( x 0 ) = Df ( y 0 ) Dg ( x 0 )
MATLAB will directly apply the chain rule when instructed to differentiate composite functions.
110
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-2
dg
Consider the functions f ( x , y ) = x 2 + y and h (u) = [sin(3u), cos(8u)] . If g (u) = f [h (u)] calculate at u = 0.
du
>> maple ('f:=(x,y) - > x ^ 2 + y ');
>> maple ('h: = u - > (sin(3*u), cos(8*u))');
>> maple ('g: = u - > (f@h) (u)');
>> pretty(sym(maple('diff(g(u),u)')))
6 sin(3 u) cos(3 u) - 8 sin(8 u)
This is the derivative of the composite function. Now, we find the value of its derivative at u = 0, defining it first as
a function of u:
>> maple ('dcomp: = u - > diff (g (u) u)');
dcomp: = u - > diff (g (u), u)
>> numeric (maple ('subs (u = 0, dcomp (u))'))
ans =
0
EXERCISE 5-3
¶z ¶z
Calculate and given that:
¶x ¶y
u2 + v2
z= , u = e x - y and v = e xy .
u2 - v2
>> maple ('f:=(u,v) - >(u^2+v^2) /(u^2-v^2)');
Now we define the following vector function:
>> maple ('g:=(x,y) - > (exp(x-y), exp(x*y))');
We define the function z as the composite of the previous two functions:
>> maple ('z:=(x,y) - > (f@g)(x,y)');
111
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
Finally, we differentiate the composite function. MATLAB automatically differentiates composite functions
via the usual syntax:
>> pretty(sym(maple('simplify(diff(z(x,y),x))')))
(- 1 + y) exp(2 x - 2 y + 2 x y)
4 --------------------------------
2
(exp(2 x - 2 y) - exp(2 x y))
>> pretty(sym(maple('simplify(diff(z(x,y),y))')))
(1 + x) exp(2 x - 2 y + 2 x y)
4 -------------------------------
2
(exp(2 x - 2 y) - exp(2 x y))
Alternatively we could have directly calculated the composite function and differentiated it as follows:
>> pretty(sym(maple('simplify(f(g(x,y)))')))
exp(2 x - 2 y) + exp(2 x y)
---------------------------
exp(2 x - 2 y) - exp(2 x y)
Now, the problem would be reduced to finding the partial derivatives of this expression with respect to x and y,
which is trivial with MATLAB.
>> maple('comp:=(x,y)->simplify(f(g(x,y)))');
>> pretty(sym(maple('simplify(diff(comp(x,y),x))')))
(- 1 + y) exp(2 x - 2 y + 2 x y)
4 --------------------------------
2
(exp(2 x - 2 y) - exp(2 x y))
>> pretty(sym(maple('simplify(diff(comp(x,y),y))')))
(1 + x) exp(2 x - 2 y + 2 x y)
4 -------------------------------
2
(exp(2 x - 2 y) - exp(2 x y))
Obviously, the result obtained is the same.
112
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
¾ [ F1 ( x , y ),..., Fm ( x , y )]
( x , y ) ¾®
F
If Fi (i = 1, 2,..., m) are differentiable with continuous derivatives up to order r and the Jacobian matrix
J = ∂ (F1,..., Fm) / ∂ (y1,..., ym) has non-zero determinant at a point ( x 0 , y 0 ) such that F ( x 0 , y 0 ) = 0 , then there is an
open U Ì Rn containing x 0 , an open V Ì Rm containing y 0 and a single-valued function f : U → V such that
F éë x , f ( x )ùû = 0 "x Î U and f is differentiable of order r with continuous derivatives.
This theorem guarantees the existence of certain derivatives of implicit functions. MATLAB allows differentiation
of implicit functions and offers the results in those cases where the hypothesis of the theorem are met.
EXERCISE 5-4
Find the conditions on x, y, z so that the surface defined by the following equation defines an implicit function
z = z(x,y) where z is a differentiable function:
x 3 + 3 y 2 + 8 xz 2 - 3 yz 3 = 1
Calculate:
¶z ¶z ¶ 2 z ¶2z
, , 2 and 2 .
¶x ¶y ¶x ¶y
>> maple('f:=(x,y,z)->x^3+3*y^2+8*x*z^2-3*z^3*y-1');
In order for the hypotheses of the implicit function theorem to be met, the partial derivative of f with respect to the
variable z must be non-zero.
>> pretty(sym(maple('diff(f(x,y,z),z)')))
2
16 x z - 9 z y
Thus the required condition is 16xz ¹ 9yz 2.
Now we can calculate the partial derivatives, assuming that the above condition is met.
The implicit (total) derivative of f and D(z) are found as follows:
>> pretty(sym(maple('D(f(x,y,z)=0)')))
2 2 2 3
3 D(x) x + 6 D(y) y + 8 D(x) z + 16 x D(z) z - 9 D(z) z y - 3 z D(y) = 0
>> pretty(sym(maple('D(z):=solve(D(f(x,y,z)=0),D(z))')))
>> pretty(sym(maple('D(z)')))
2 2 3
3 D(x) x + 6 D(y) y + 8 D(x) z - 3 z D(y)
D(z) := - ---------------------------------------------
2
16 x z - 9 z y
113
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
As y does not depend on x (x and y are, by hypothesis, independent variables, and the variable z depends on
x and y ), we put D(y) = 0 in the expression for D(z)/D(x):
>> pretty(sym(maple('simplify(subs(D(y)=0,D(z)/D(x)))')))
2 2
3 x + 8 z
- -------------------
z (16 x - 9 z y)
This gives the expression for ∂z/∂x. For ∂z/∂y, we similarly calculate the following:
>> pretty(sym(maple('simplify(subs(D(x)=0,D(z)/D(y)))')))
3
2 y - z
- 3 ------------------
z (16 x - 9 z y)
To calculate ∂ 2z /∂x 2, we find the partial derivative of ∂z/∂x with respect to x:
>> pretty(simple(sym(diff(maple('simplify(subs(D(y)=0,D(z)/D(x)))'),'x'))))
2 2
24 x - 27 x z y - 64 z
- 2 ----------------------------
2
z (16 x - 9 z y)
To calculate ∂ 2z /∂y 2, we find the partial derivative of ∂z/∂y with respect to y:
>> pretty(simple(sym(diff(maple('simplify(subs(D(x)=0,D(z)/D(y)))'),'y'))))
4
32 x - 9 z
- 3 -------------------
2
z (16 x - 9 z y)
EXERCISE 5-5
Show that near the point (x, y, u, v ) = (1,1,1,1) the system
xy + yvu 2 = 2
xu 3 + y 2v 4 = 2
114
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
The functions are differentiable and have continuous derivatives. We need to show that the corresponding
Jacobian determinant is non-zero at the point (1,1,1,1):
>> maple('f1:=(x,y,u,v)-> x*y +y*v*u^2-2');
>> maple('f2:=(x,y,u,v)->x*u^3+y^2*v^4-2');
>> sym(maple('jacobian([f1(x,y,u,v),f2(x,y,u,v)],[u,v])'))
ans =
[2*y*v*u, y*u^2]
[3*x*u^2, 4*y^2*v^3]
>> J=sym(maple('subs(x=1,y=1,u=1,v=1,jacobian([f1(x,y,u,v),f2(x,y,u,v)],[u,v]))'))
J =
[2, 1]
[3, 4]
>> determ(J)
ans =
5
Thus the assumptions of the implicit function theorem are met and the proposed system can be solved uniquely.
115
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-6
116
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
Observe that this is indeed the reciprocal of the determinant found above.
Next we find the value of the inverse Jacobian at the point (p/4,−p/4) :
>> numeric(subs(subs(determ(I),pi/4,'x'),-pi/4,'y'))
ans =
0.3821
>> numeric(subs(subs(symdiv(1,determ(J)),pi/4,'x'),-pi/4,'y'))
ans =
0.3821
Here it is clear that the determinant of the Jacobian of the inverse function is the reciprocal of the determinant of
the Jacobian of the function.
EXERCISE 5-7
Demonstrate that the transformation between cartesian and polar coordinates complies with the assumptions of
the inverse function theorem.
We know that the transformation equations are:
x = a cos[b], y = a sin[b]
Obviously, the functions have continuous partial derivatives. Let's see if the determinant of the Jacobian of the
transformation is non-zero:
>> J = simple (sym (maple ('jacobian ([a * cos (b), a * sin (b)], [a, b])')))
J =
[cos (b),-a * sin (b)]
[sin(b), a*cos(b)]
>> pretty(simple(sym(determ(J))))
a
We see that the Jacobian of the transformation is non-zero (a ¹ 0). Thus, the inverse function theorem is
applicable. The determinant of the Jacobian of the inverse transformation will be 1/a.
117
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-8
Consider the function f ( x , y ) = e x - y and the transformation u = u(x,y ) = x + y, v = v(x,y ) = x. Calculate f (u,v ).
>> maple ('f:=(x,y) - > exp(x-y)');
This transformation fulfils the conditions of the inverse function theorem (continuous partial derivatives and a
Jacobian with non-zero determinant).
We find the inverse transformation and its Jacobian to apply the change of variables theorem:
>> maple ('solve({u=x+y,v=x},{x,y})')
ans =
{x = v, y = u-v}
>> sym (maple ('jacobian([v,u-v],[u,v])'))
ans =
[0, 1]
[1, - 1]
>> pretty (simple (sym (maple ('f (v, u-v) * abs (det (jacobian([v,u-v],[u,v])))'))))
exp(2 v-u)
The requested function is f(u,v)= e2v-u.
118
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
Here
x = ( x1 , x 2 x n ), a = (a1 , a2 an ),ti = xi - ai (i = 1, 2 ,n)
R = remainder.
Normally, the series are given up to order 2.
EXERCISE 5-9
Find the Taylor series of order 2 at the point (1,0) of the function:
2
f ( x , y ) = e ( x -1) cos( y ).
>> maple('f:=(x,y)->exp((x-1)^2)*cos(y)');
>> pretty(sym(maple('f(1,0)+D[1](f)(1,0)*(x-1)+D[2](f)(1,0)*(y)+
+(1/2!)*(D[1,1](f)(1,0)*(x-1)^2+ 2*D[1,2](f)(1,0)*( x-1)*y+D[2,2](f)(1,0)*(y)^2)')))
2 2
1 + (x - 1) - 1/2 y
EXERCISE 5-10
119
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
Definition of a scalar potential of a vector field: A vector field F is called conservative if there is a differentiable
function f such that F = Df . The function f is known as a scalar potential function for F .
Definition of the curl of a vector field: The curl of a vector field F ( x , y , z ) = Mi + Nj + Pk is the following:
æ ¶P ¶N ö æ ¶P ¶M ö æ ¶N ¶M ö
curl F ( x , y , z ) = D ´ F ( x , y , z ) = ç - ÷i - ç - ÷ j +ç - ÷k
è ¶y ¶z ø è ¶x ¶z ø è ¶x ¶y ø
Definition of a vector potential of a vector field: A vector field F is a vector potential of another vector field
G if F = curl (G).
Definition of the divergence of a vector field: The divergence of the vector field F ( x , y , z ) = Mi + Nj + Pk is the
following:
¶M ¶N ¶P
diverge F ( x , y , z ) = D × F ( x , y , z ) = + +
¶x ¶y ¶z
Definition of the Laplacian: The Laplacian is the differential operator defined by:
¶2 ¶2 ¶2
Laplacian = D 2 = D × D = + +
¶x 2 ¶y 2 ¶z 2
120
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-11
x y z
[---------------------, ----------------------, ---------------------]
2 2 2 3/2 2 2 2 3/2 2 2 2 3/2
(1 - x - y - z ) (1 - x - y - z ) (1 - x - y - z )
>> pretty(simple(sym(maple('laplacian(w,[x,y,z])'))))
3
---------------------
2 2 2 5/2
(1 - x - y - z )
EXERCISE 5-12
ans =
true
The response true is given because there is a scalar potential, i.e., because the curl of the given vector
field is null:
>> sym (maple ('curl([1,2*y,3*z^2],[x,y,z])'))
ans =
[0, 0, 0]
The value of the scalar potential is contained in the variable P :
>> pretty(sym(maple('P')))
2 3
y + x + z
121
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
If we calculate the gradient of this scalar potential, we get the initial vector field.
>> pretty(simple(maple('grad(y^2+x+z^3,[x,y,z])')))
2
[1, 2 y, 3 z ]
EXERCISE 5-13
EXERCISE 5-14
ans =
true
Thus, we know that there is a vector potential, and the divergence of the vector field is zero:
>> sym (maple ('diverge([x*z,-y*z,y],[x,y,z])'))
ans =
0
122
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
2 2
[- 1/2 y (z + y), - 1/2 x z , 0]
If we calculate the curl of the previous result, we obtain the initial vector field:
>> pretty (sym (maple ('curl ([-1/2 * y *(z^2+y), - 1/2 * x * z ^ 2, 0], [x, y, z])')))
[x z,- y z, y]
Now, we calculate the vector potential of the other two fields:
>> pretty(simple(sym(maple('print(Q)'))))
2 2
[- 1/2 y (z + y), - 1/2 x z , 0]
[-z sin (y) - y sin (x), - x z cos(y), 0]
If we calculate the curl of these vector potentials, we obtain the initial vector fields:
>> pretty(sym(maple('curl([- sin(z)*x, - x*exp(y)*z, 0],[x,y,z])')))
123
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
Rectangular to cylindrical:
r = x2 + y2
y
q = tan -1
x
z=z
Spherical to rectangular:
x = r sin f cosq
y = r sin f sin q
z = r cosf
Rectangular to spherical:
r = x2 + y2 + z2
y
q = tan -1
x
z
f = cos -1
x + y2 + z2
2
EXERCISE 5-15
Express the point given in spherical coordinates by (4,p/6,p/4) in rectangular and cylindrical coordinates.
>> maple ('sphrec(4,pi/6,pi/4)')
ans =
[2 ^(1/2) * 3 ^(1/2), 2 ^(1/2), 2 * 2 ^(1/2)]
124
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-16
Express the point given in cylindrical coordinates by (4,5p/6,3) in rectangular and spherical coordinates.
>> maple('cylrec(4,5*Pi/6,3)')
ans =
[- 2 * 3 ^ (1/2), 2, 3]
Now, we convert the rectangular to spherical coordinates:
>> maple ('recsph (- 2 * 3 ^ (1/2), 2, 3)')
ans =
5, - 1/6 * pi, acos(3/5)
EXERCISE 5-17
Express the equations of the surfaces xz = 1 and x 2 + y 2 + z 2 = 1 in spherical coordinates.
>> pretty(sym(maple('subs(x=r*sin(a)*cos(t),y=r*sin(a)*sin(t),z=r*cos(a), x*z=1)')))
2
r sin(a) cos(t) cos(a) = 1
>> pretty(simple(sym(maple('subs(x=r*sin(a)*cos(t),y=r*sin(a)*sin(t),
z=r*cos(a), x^2+y^2-z^2=1)'))))
2
- r cos(2 a) = 1
125
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
EXERCISE 5-18
Express in cartesian and spherical coordinates the equation of the surface which in cylindrical coordinates is
given by z = r 2 (1 + sin (t )).
>> pretty(sym(maple('expand(simplify(subs(r=sqrt(x^2+y^2),t=arctan(y/x),
z=z, z=r^2*(1+sin(t))),symbolic))')))
2 2 2 2 1/2
z = x + y + (x + y ) y
Now let us convert the cartesian to spherical coordinates:
>> pretty(sym(maple('assume(t,real):expand(simplify(subs(x=r*sin(a)*cos(t),
y=r*sin(a)*sin(t), z=r*cos(a), z=(x^2+y^2)*(1+sin(atan(y/x)))),symbolic))')))
2 2 2 2 2 2
r cos(a) = r sin(t) - r cos(a) sin(t) + r - r cos(a)
EXERCISE 5-19
Find the unit tangent, the unit normal, and the unit binormal vectors of the twisted cubic: x = t, y = t 2, z = t 3.
We begin by restricting the variable t to the real field:
>> maple('assume(t,real)') ;
We define the symbolic vector V as follows:
>> syms t, V=[t,t^2,t^3]
V =
[t, t ^ 2, t ^ 3]
The tangent vector is calculated by:
>> tang= diff (V)
tang =
[1, 2 *, 3 * t ^ 2]
126
www.pdfgrip.com
Chapter 5 ■ Vector Differential Calculus and Theorems in Several Variables
un =
[ (-2*t-9*t^3)/(9*t^4+9*t^2+1)^(1/2)/(1+4*t^2+9*t^4)^(1/2),
(1-9*t^4)/(9*t^4+9*t^2+1)^(1/2)/(1+4*t^2+9*t^4)^(1/2),
(6*t^3+3*t)/(9*t^4+9*t^2+1)^(1/2)/(1+4*t^2+9*t^4)^(1/2)]
The unit binormal vector is the vector product of the tangent vector and the unit normal vector.
>> ub=simple(cross(ut,un))
ub =
[3*t^2/(9*t^4+9*t^2+1)^(1/2),-3*t/(9*t^4+9*t^2+1)^(1/2),1/(9*t^4+9*t^2+1)^(1/2)]
The unit binormal vector can also be calculated via (v'∧v ") / |v'∧v" | as follows:
>> ub=simple(v1/sqrt(dot(v1,v1)))
ub =
[3*t^2/(9*t^4+9*t^2+1)^(1/2),-3*t/(9*t^4+9*t^2+1)^(1/2),1/(9*t^4+9*t^2+1)^(1/2)]
We have calculated the Frenet frame for a twisted cubic.
127
www.pdfgrip.com
Chapter 6
MATLAB works with integral calculus in a clear and simple way. The number of functions that enables you to work in
this area is not very high, but they are very efficient in solving integration problems. You can calculate the indefinite
integral of most integrable functions whose structure is not very complicated; for example, functions that involve
simple logarithms, exponentials, rational, trigonometric functions, inverse trigonometric functions, etc. Definite and
improper integrals do not present problems for MATLAB. Double integrals, triple integrals and n-fold integrals are
also easily found.
The most commonly used MATLAB commands for integral calculus are:
syms x, int(f(x), x) or int(‘f(x)’, ‘x’)
Computes the indefinite integral ò f ( x )dx
int (int (‘f(x,y)’, ‘x’), ‘y’)
Computes the double integral ò ò f ( x , y )dxdy
syms x y int (int (f(x,y), x), y)
Computes the double integral òò f ( x , y )dxdy
int(int(int(...int(‘f(x,y...z)’, ‘x’), ‘y’)...), ‘z’)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz
syms x y z, int(int(int(...int(f(x,y...z), x), y)...), z)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz
syms x a b, int(f(x), x, a, b)
b
Computes the definite integral òa
f ( x )dx
int (‘f (x)’, ‘x’, ‘a’, ‘b’)
b
Computes the definite integral òa
f ( x )dx
int(int(‘f(x,y)’, ‘x’, ‘a’, ‘b’), ‘y’, ‘c’, ‘d’))
b d
Computes the definite double integral ò ò
a c
f ( x , y )dxdy
syms x y b c d, int (int (f(x,y), x, a, b), y, c, d)
b d
Computes the definite double integral ò ò
a c
f ( x , y )dxdy
int(int(int(....int(‘f(x,y,...,z)’, ‘x’, ‘a’, ‘b’), ‘y’, ‘c’, ‘d’),...), ‘z’, ‘e’, ‘f’)
b d f
Computes ò ò ¼ ò f ( x , y ,¼ , z ) dxdy ¼ dz
a c e
129
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
In addition, with the prior use of the command maple, the following commands can be used for integration:
int (f (x), x) Computes the indefinite integral ò f ( x )dx
b
int (f (x), x = a...b) Computes the definite integral òa
f ( x )dx
int(f(x), x = a..b, continuous)
Computes the definite integral removing the check for continuity of the
function in the interval (a, b)
residue (f(x), x = a)
Returns the algebraic residue of f(x) at the point a. Used to solve
integrals by the method of residues. The algebraic residue is the
coefficient of x -1 in the Laurent expansion of f(x) at a.
readlib (singular):singular(f(x))
Returns the singularities of the function f(x)
readlib (singular): singular (f (x), x)
Returns the singularities of the function f(x) treating all variables other
than x as constants
readlib (singular): singular (f(x, y, z,...), {x, y,z,...})
Returns the singularities of the function f(x,y,z,...) treating all variables
other than x, y,z,... as constants
op (Int (f (x), x))
Shows the internal structure of the Int command and allows access to
its operands
op (n, Int (f (x), x)
Removes the n-th operand of int. The first operand is the integrand, the
second operand is the variable of integration or x = a...b for definite
integrals, and the third operand is continuous, if it is being used. The
operands can be replaced with the command subsop.
int (int (f(x,y), x), y) computes the indefinite integral ò ò f ( x , y )dxdy
int(int(int(...int(f(x,y,...,z), x), y),...), z)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz
int(int(f(x,y),x=a..b),y=c..d)
b d
Computes the definite integral ò ò
a c
f ( x , y )dxdy
int (int (int (... int (f(x,y,...,z), x = a...b), y = c...d),...), z = e...f)
b d f
Computes ò ò ¼ ò f ( x , y ,¼ , z ) dxdy ¼ dz
a c e
130
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
(n + 1)
x
---------
n + 1
We have seen how MATLAB allows you to include parameters in the integrals, which are treated as generic
constants:
>> pretty(simple(int(int('x^3+y^3','x',0,'a'),'y',0,'b')))
4 4
1/4 a b + 1/4 a b
>> syms x y a b, pretty(simple(int(int(x^3+y^3,x,0,a),y,0,b)))
4 4
1/4 a b + 1/4 a b
>> pretty(simple(int(int('x^3+y^3','x',0,'a'),'y',0,'x')))
4 4
1/4 a x + 1/4 a x
In these last examples definite integrals have been solved, one of which has a variable limit of integration,
which is perfectly valid with MATLAB.
131
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-1
>> pretty(simple(int('sec(x)*csc(x)')))
log(tan(x))
>> pretty(simple(int('x*cos(x)')))
cos(x) + x sin(x)
>> pretty(simple(int('acos(2*x)')))
2 1/2
x acos(2 x) - 1/2 (1 - 4 x )
Exercise 6-2
132
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-3
æxö
tan -1 ç ÷
è2ø a sin( x ) cos3 ( x )
ò x 2 + 4 dx , ò (1 - x 2 )3/2 dx , ò sin( x ) dx.
>> maple('with(student)');
>> pretty(simple(sym(maple('changevar(t=atan(x/2),int(atan(x/2)/(4+x^2),x),t)'))))
2 2
- 1/16 log(1 + 1/2 i x) + 1/8 log(2) log(1/16 x + 1/4)
2
- 1/8 dilog(1/2 + 1/4 i x) - 1/16 log(1 - 1/2 i x)
- 1/8 dilog(1/2 - 1/4 i x)
133
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Results can easily be found by using the integral operator in inert mode, i.e. Int, making the substitution and
finding the value of the integral later with value as follows:
>> pretty(sym(maple('value(simplify(changevar(t=arctan(x/2),
Int(arctan(x/2)/(4+x^2),x),t)))')
2
1/4 t
Now we undo the change of variable to obtain the final result:
>> pretty(subs('1/4*t^2','atan(x/2)'))
2
1/4 atan(1/2 x)
>> pretty(simple(sym(maple('changevar(t=asin(x), Int(asin(x)/(1-x^2)^(3/2),x),t)'))))
t tan(t) + log(cos(t))
Now we reverse the change of variable:
>> pretty(subs('t*tan(t)+t*log(t)','asin(x)'))
asin(x) x
----------- + asin(x) log(asin(x))
2 1/2
(1 - x )
>> pretty(simple(sym(maple('changevar(t=sin(x), Int(cos(x)^3/sin(x)^(1/2),x),t)'))))
1/2 2
- 2/5 t
- 5) (t
>> pretty(subs('-2/5*t^(1/2)*(t^2-5)','sin(x)'))
1/2 2
- 2/5 sin(x) (sin(x) - 5)
134
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-4
3
- 1/3 atanh(-------------)
2 1/2
(9 - 4 x )
For the third integral, we use the change of variable 3 + 5x 3 = t 4:
>> pretty(simple(sym(maple('changevar(t=(3+5*x^3)^(1/4),int(x^8* (3+5*x^3)^(1/4),x),t)'))))
3 6 9 3 1/4
4/73125 (288 - 120 x + 125 x + 1875 x ) (3 + 5 x )
135
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
MATLAB incorporates the command maple(‘intparts’) which allows you to find an integral directly by parts.
This command belongs to the student library, which has to be previously loaded into the memory with the command
maple (‘with(student)’). The syntax is as follows:
maple (‘intparts (int (f (x), x), u (x))’)
Finds the integral of f(x) by parts, where u(x) is a factor of f(x) which
will be differentiated when calculating the integral.
As an example, we calculate the integral of the function (3x2 + 2x − 7)cos (x):
>> maple('with(student)');
>> pretty(simple(sym(maple('intparts(int((3*x^2+2*x-7)*cos(x),x),3*x^2+2*x-7)'))))
2
3 x sin(x) - 13 sin(x) + 6 cos(x) x + 2 cos(x) + 2 x sin(x)
Exercise 6-5
136
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-6
ò sin
13
( x )cos15 ( x ) dx , òx 7 (m + nx )1/2 dx , òe 12 x cos11 ( x ) dx .
137
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
For a planar curve with parametric coordinates x = x (t) and y = y (t), the arc length of the curve between the
points corresponding to the values t = t0 and t = t1 of the parameter is given by the expression:
t1
L=ò x ¢(t )2 + y ¢(t )2 dt
t0
For a curve in polar coordinates with equation r = f (a), the arc length of the curve between the points
corresponding to the parameter values a = t0 and a = t1 is given by the expression:
a1
L=ò r 2 + r ¢(a )2 dr
a0
For a space curve with parametric coordinates x = x (t), y = y (t), z = z (t), the arc length of the curve between the
points corresponding to the parameter values t = t0 and t = t1 is given by the expression:
t1
L=ò x ¢(t )2 + y ¢(t )2 + z ¢(t )2 dt
t0
For a space curve in cylindrical coordinates with equations x = r · cos(a), y = r · sin (a), z = z, the arc length of the
curve between the points corresponding to the parameter values a = a0 and a = a1 is given by the expression:
a1
L=ò r 2 + r ¢2 + z ¢2 dr
a0
138
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
For a space curve in spherical coordinates with equations given by x = r · sin(a) · cos(b), y = r · sin(a) · sin(b),
z = r · cos(a), the arc length of the curve between the points corresponding to the values a = a0 and a = a1 is given by
the expression:
a1
L=ò dr 2 + r 2 da 2 + r 2 sin 2 (a )db 2 dr 2
a0
Exercise 6-7
An electrical cable hangs between two towers which are 80 meters apart. The cable adopts the shape of a
catenary whose equation is:
x
y = 100 cosh
100
Calculate the arc length of the cable between the two towers.
Given the symmetry of the curve with respect to the y-axis, the limits of integration will be –40 and 40:
>> f ='100*cosh(x/100)';
>> diff(f,'x')
ans =
sinh(1/100*x)
>> pretty(simple(int('(1+sinh(1/100*x)^2)^(1/2)','x','-40','40')))
100 exp (2/5) - 100 exp(-2/5)
If we want the approximate result, we use the command numeric:
>> numeric(simple(int('(1+sinh(1/100*x)^2)^(1/2)','x','-40','40')))
ans =
82.1505
Result: 82.1505 m.
139
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-8
Calculate the arc length of the space curve represented by the parametric equations:
2
4 t
x = t y = t 3/2 z = 6
3 2
between t = 0 and t = 2.
>> pretty(simple(int('(diff(t,t)^2+diff((4/3)*t^(3/2),t)^2+diff((1/2)*t^2,t)^2)^(1/2)',
't','0','2')))
1/2 1/2
2 13 - 3/2 log(13 + 4) - 1 + 3/2 log(3)
Now we approximate the result:
>> numeric(simple(int('(diff(t,t)^2+diff((4/3)*t^(3/2),t)^2+diff((1/2)*t^2,t)^2)^(1/2)',
't','0','2')))
ans =
4.8157
Result: 4.8157.
Exercise 6-9
Find the arc length for values from a = 0 to a = 2pi of the cardioid given in polar coordinates by r = 3-3cos(a).
>> r = '3-3*cos(a)';
>> diff(r,'a')
ans =
3 * sin (a)
>> R = simple (int ('((3-3 * cos (a)) ^ 2 + (3 * sin (a)) ^ 2) ^(1/2) ',' a ', ' 0','2 * pi'))
R=
24
Result: 24.
140
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
where x = a and x = b are the abscissas of the end points of the curve.
If the curve is given in parametric coordinates x = x (t) and y = y(t), the area is given by the integral:
b
S= òa
y (t ) x ¢(t ) dt
where values of the parameter t = a and t = b correspond to the end points of the curve.
If the curve is given in polar coordinates r = f (a), the area is given by the integral:
1 a1
2 òa0
S= f (a )2 da
for the parameter values a = a0 and a = a1 corresponding to the end points of the curve.
To calculate the area between two curves with equations y = f(x) and y = g(x), we use the integral:
b
S = ò f ( x ) - g ( x ) dx
a
where x = a and x = b are the abscissas of the end points of the two curves.
When calculating these areas it is very important to take into account the sign of the functions involved since the
integral of a negative portion of a curve will be negative. One must divide the region of integration so that positive and
negative values are not computed simultaneously. For the negative parts one takes the modulus.
Exercise 6-10
Find the area of the region bounded by the curves defined below:
f ( x ) = 2 - x 2 and g ( x ) = x .
We will have to find the abscissas of the points where the curves meet. Consider the equation:
2 - x2 = x
>> solve('2-x^2=x')
ans =
[-2]
[1]
Now, we can already find the area that we have to calculate:
>> int('2-x^2-x','x',-2,1)
ans =
9/2
141
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
>> fplot('[2-x^2,x]',[-2,1])
Figure 6-1.
Exercise 6-11
Calculate the area under the normal curve between the limits -1.96 and 1.96.
2
e x /2
1.96
We need to calculate the integral ò-1.96 2p dx .
numeric(int('exp(-x^2/2)/(2*pi)^(1/2)','x',-1.96,1.96))
ans =
0.9500
142
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-12
Calculate the area of the ellipse with major and minor axes of length a and b.
2 2
y
As is well known, the equation of this ellipse is: x 2 + 2 = 1
a b b 2
Rearranging this equation for y and assuming the positive root, we obtain y = a - x2 .
a
By the symmetry of the ellipse, the area will be four times the integral of this expression between 0 and a.
a 4b a 2 - x 2
Therefore we calculate: ò
0 a
dx .
>> pretty(int('(4*b*(a^2-x^2)^(1/2))/a','x',0,'a'))
pi a b
The requested result is therefore pab.
Exercise 6-13
Calculate the length of the perimeter and the area enclosed by the curves:
cos ( t )[2 - cos(2t )] sin ( t )[2 + cos(2t )]
x (t ) = , y (t ) = .
4 4
First, in order to get an idea of what the problem is, we graphically represent the curve (see Figure 6-2):
>> x = (0:.1:2*pi);
>> t = (0:.1:2*pi);
>> x = cos(t).*(2-cos(2*t))./4;
>> y = sin(t).*(2+cos(2*t))./4;
>> plot(x,y)
143
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Figure 6-2.
We see that as the parameter ranges between 0 and 2p the curve described is closed. We can then calculate its
length and the area that it encloses.
>> A = simple (diff ('cos (t) * (2-cos(2*t)) / 4'))
A =
-3/8 * sin (t) + 3/8 * sin(3*t)
>> B = simple (diff (' sin (t) * (2 + cos(2*t)) / 4'))
B =
3/8 * cos (t) + 3/8 * cos(3*t)
Now we apply the formula for the length of a curve in parametric coordinates (considering only half a quadrant):
>> int ('sqrt ((-3/8 * sin (t) + 3/8 * sin(3*t)) ^ 2 + (3/8 * cos (t) + 3/8 * cos(3*t)) ^ 2)',
't', 0, pi/4)
ans =
3/8
144
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-14
Calculate the length and the area enclosed by each of the following curves given in polar coordinates:
r = cos(2a ) and r = sin ( 2a ) .
145
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Figure 6-3.
Now, we calculate the arc lengths and the areas enclosed by both curves. The first curve repeats the structure for
a between 0 and p/4 four times, and the second for a between 0 and p/2:
>> r = 'sqrt(cos(2*a))';
>> A = simple(diff(r,'a'))
A =
-1/cos(2*a) ^(1/2) * sin(2*a)
We calculate the area S1 enclosed by the first curve. The area under the curve in the first quadrant is:
>> I = simple(int('(sqrt(cos(2*a)))^2','a',0,pi/4))
I=
1/2
146
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
147
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-15
Calculate the area enclosed by the curve y = sin(x) and y = cos(x) for x varying between 0 and 2p.
Figure 6-4 shows the graph of both curves on the same axis, so we can see the points of intersection and the
region between the curves.
>> fplot ('[sin (x), cos (x)], [0, 2 * pi])
Figure 6-4.
We need to know the coordinates of the points of intersection of the two curves. To do this, we solve the system
formed by their two defining equations. The command solve does not solve the problem. Therefore, since the
graph indicates that there are solutions in the intervals (0,1) and (1,2p), we will use the command maple(‘fsolve’)
to find these solutions:
>> maple ('fsolve (sin (x) = cos (x), x, 0.. 1)')
ans =
.7853981633974483
>> maple ('fsolve (sin (x) = cos (x), x, 1.. 2 * pi)')
ans =
3.926990816987242
148
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
where x = a and x = b are the x-coordinates of the end points of the revolving curve. If the curve is given in parametric
coordinates x = x(t) and y = y(t) and the end points of the curve correspond to t = t0 and t = t1 then the rotational surface
area is given by the integral:
t1
S = 2p ò y (t ) x ’(t ) + y ’( x )2 dt
t0
The area of the surface generated by rotating the curve with equation x = f (y) around the y axis is given by
the integral:
b
S = 2p ò f ( y ) 1 + f ’( y )2 dy
a
where y = a and y = b are the y-coordinates of the end points of the rotating curve.
149
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-16
Calculate the area of the surface generated by rotating the cubic y = 12x − 9x 2 + 2x 3 around the OX axis, for x
between 0 and 5/2.
The surface of revolution has equation y 2 + z 2 = (12x − 9x 2 + 2x 3)2, and for the purposes of producing the graph
shown in Figure 6-5, we can parameterize the surface as follows:
150
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-17
Calculate the area of the surface of revolution given by rotating the parametric curve defined by x (t ) = t − sin(t )
and y (t ) = 1 − cos(t ) around the OX axis, for t between 0 and 2p .
>> x = 't - sin (t)';
>> y = '1 - cos (t)';
>> [diff(x),diff(y)]
ans =
[1-cos (t), sin (t)]
Now, we find the requested surface area by integration:
>> I = numeric(int('(1-cos(t))*sqrt(abs((1-cos(t))^2+sin(t)^2))','t',0,2*pi))
I =
10.6667
The area will be 2p(10.6667) = 21.3334p square units.
where x = a and x = b are the x-coordinates of the end points of the rotating curve.
The volume generated by rotating the curve with equation x = f (y) around the y-axis is given by the integral:
b
V = p ò f ( y )2 dy
a
where y = a and y = b are the y-coordinates of the end points of the rotating curve.
If one cuts a volume by planes parallel to one of the three coordinate planes (for example, the plane z = 0) and if
the equation S(z) of the curve given by the cross section is given in terms of the distance of the plane from the origin
(in this case, z) then the volume is given by:
z2
V = ò S (z ) dz
z1
151
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-18
x2 y2
Calculate the volume generated by rotating the ellipse + = 1 around the OX axis and around the OY axis.
4 9
We depict half of each of the generated surfaces of revolution next to each other in Figure 6-6. The final volume
will therefore be twice that of the represented figure. The equation of the surface of revolution around the
æ x2 ö
OX axis is y 2 + z 2 = 9 ç 1 - ÷ , which can be parameterized as follows:
è 4 ø
1 1
æ t 2 ö2 æ t 2 ö2
x = t , y = 3 cos ( u ) ç 1 - ÷ , z = 3 sin(u) ç 1 - ÷
è 4ø è 4ø
Figure 6-6.
æ y2 ö
The equation of the surface of revolution around the OY axis is x 2 + z 2 = 4 ç 1 - ÷ , which can be parameterized
è 9 ø
as follows:
1 1
æ t 2 ö2 æ t 2 ö2
x = 2 cos ( u ) ç 1 - ÷ , y = t , z = 2 sin(u) ç 1 - ÷
è 9ø è 9ø
>> t = (0:.1:2);
>> u = (0:.5:2*pi);
>> x = ones (size (u))'* t;
>> y = cos(u)'*3*(1-t.^2/4).^(1/2);
>> z = sin(u)'*3*(1-t.^2/4).^(1/2);
>> subplot(1,2,1)
>> surf(x,y,z)
152
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
>> subplot(1,2,2)
>> x = cos(u)'*2*(1-t.^2/9).^(1/2);
>> y = ones (size (u))'* t;
>> z = sin(u)'*2*(1-t.^2/9).^(1/2);
>> surf (x, y, z)
Now we calculate the generated volumes by integration:
>> V1 = int('pi*9*(1-x^2/4)','x',-2,2)
V1 =
24*pi
>> V2 = int('pi*4*(1-y^2/9)','y',-3,3)
V2 =
16*pi
Exercise 6-19
Calculate the volume generated by rotating the curve given in polar coordinates by r = 1 + cos(a) about
the OX axis.
The curve is given in polar coordinates, but there is no problem in calculating the values of x(a) and y(a) needed to
implement the volume formula in cartesian coordinates. We have:
x (a) = r (a) cos (a) = (1 + cos (a)) cos (a)
y (a) = r (a) sin (a) = (1 + cos (a)) sin (a)
>> x ='(1 + cos (a)) * cos (a)';
y ='(1 + cos (a)) * sin (a)';
>> A = simple(diff(x))
A =
-sin(2*a) - sin (a)
>> I = int ('((1 + cos (a)) * sin (a)) ^ 2 * (- sin(2*a) - sin (a))', 'a', 0, pi)
I =
-8/3
The required volume is p (abs (I)) = 8p/3 cubic units.
153
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
ò F • ds = ò F [c(t )] • c¢(t ) dt
b
c a
Exercise 6-20
Consider the curve c (t ) = [sin (t ) ,cos (t ) ,t ] with 0 < t < 2p , and the vector field F ( x , y , z ) = xi + y j + zk .
Calculate òc F × ds.
derivC = [diff ('sin (t)'), diff('cos (t)'), diff ('t')]
derivC =
[cos (t),-sin(t), 1]
>> pretty (int ('dotprod (vector ([sin (t), cos (t), t]), vector ([cos (t),-sin(t), 1]))
', ', 0, 2 * pi))
2
2 pi
Exercise 6-21
derivC = [diff ('cos (a) ^ 3'), diff (' sin (a) ^ 3'), diff('a')]
derivC =
[- 3 * cos (a) ^ 2 * sin (a) 3 * sin (a) ^ 2 * cos (a) 1]
>> pretty (sym (numeric (int ('dotprod (vector ([sin (a), cos (a) - sin (a) ^ 3 *
cos (a) ^ 3]), vector ([- 3 * cos (a) ^ 2 * sin (a), 3 * sin (a) ^ 2 * cos (a), 1]))
^(1/3)', 'a', 0, 7 * pi/2))))
1/2
5/2 + 3i
154
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-22
Let the curve c be the circle of radius a: x 2 + y 2 = a 2. Calculate the value of ò
c
x 3dy - y 3dx .
First, the circle can be parameterized (x = a cos(t ), y = a sin(t )) and we can then calculate the integral as follows:
>> pretty(int('a*cos(t)^3*diff(a*sin(t),t)-a*sin(t)^3*diff(a*cos(t),t)','t',0,2*pi))
2
3/2 a pi
155
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
2 æ 2 æ 1 ö æ 1 öö 1 æ 1
k
1 k ö
3 ln(3) + ç å ç +i ÷ ln ç + i ÷ ÷ + ç å (1+i ) ln (1 + i ) ÷
k
6 3 çè i =1 è 2 ø è 2 ø ÷ø 3 è i =1 ø
>> maple('trapezoid(x^k*ln(x), x = 1..3)')
ans =
1/2*Sum((1+1/2*i)^k*log(1+1/2*i),i = 1 .. 3)+1/4*3^k*log(3)
This expression can be written as:
1 æ 3 æ 1 ö æ 1 öö 1 k
k
çç å ç 1 + i ÷ ln ç 1 + i ÷ ÷÷ + 3 ln(3)
2 è i =1 è 2 ø è 2 ø ø 4
Exercise 6-23
The first of the integrals presents a singularity at x = 0, so there may be problems in a neighborhood of 0.
We have:
>> pretty(sym(limit(int(1/sqrt(x),x,a,b),a,0)))
1/2
2 b
>> pretty(simple(sym(int(1/sqrt(x),x,0,b))))
1/2
2 b
156
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
i.e. integration under the integral sign is valid, and the order in which the integrals are evaluated can be reversed
without affecting the result.
157
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-24
ò0 x (1 + x 2 )
dx , ò0 x 2 dx.
For the first integral, we will start by integrating the derivative of the integrand with respect to the parameter a,
which will be easier to integrate. Once this integral is found we integrate with respect to a to find the original
integral.
>> maple('assume(a>0)');
>> pretty(simple(sym((int(diff(atan(a*x)/(x*(1+x^2)),a),x,0,inf)))))
pi
--------
2 a + 2
Now we integrate this function with respect to the variable a, to find the original integral:
>> pretty(simple(sym(int(pi/(2*a+2),a))))
1/2 log(2 a + 2) pi
To solve the second integral, we consider the following, using a as a parameter:
¥ 2
ó 1 - e - ax
f (a ) = ôô dx
õ0 x2
As in the first integral, we differentiate the integrand with respect to the parameter a, find the integral, and then
integrate with respect to a. The desired integral is then given by setting a = 1.
>> Ia = simple (sym (int (diff ((1-exp(-a*x^2)) / x ^ 2, a), x, 0, inf)))
Ia =
1/2/a ^(1/2) * pi ^(1/2)
>> s = simple(sym(int(Ia,a)))
s =
a ^(1/2) * pi ^(1/2)
Putting a = 1, we obtain the value p for the original integral.
158
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
159
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Exercise 6-25
To approximate the integrals, we partition the range of integration into n intervals, find the middlesum, and then
find the limit of these sums as n tends to infinity:
>> maple('with(student)')
>> pretty(sym(maple('evalf(limit(middlesum(cos(sin(x)),x=0..1,n),n=infinity))')))
0.8687400396
> evalf(limit(middlesum(sin(x)^2,x=2..5,n),n=infinity));
1.446804
> evalf(limit(middlesum(1/log(x),x=2..7,n),n=infinity));
3.711887986
> evalf(limit(middlesum(sin(x)/x,x=1..5,n),n=infinity));
0.6038481746
Alternative methods can be used. For example, to approximate the first integral, we can calculate the limit as the
number of rectangles n tends to infinity of both the upper and lower sums. If these two limits match, then the
function is integrable and the common value is the integral.
>> maple ('f: = x - > cos (sin (x))');
>> numeric(maple('limit((1-0)/n*sum(abs(f(0+(1-0)*k/n)),k=1..n),n=infinity)'))
ans =
0.8687
>> numeric(maple('limit((1-0)/n*sum(abs(f(0+(1-0)*k/n)),k=0..n-1),n=infinity)'))
ans =
0.8687
If we calculate the integral directly, we see that the values coincide:
>> numeric (int ('cos (sin (x))', 'x', 0, 1))
ans =
0.8687
160
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
Figure 6-7.
161
www.pdfgrip.com
Chapter 6 ■ Integration and Applications
162
www.pdfgrip.com
Chapter 7
In this chapter we will describe how to solve multivariate integrals with MATLAB, and give applications of double
and triple integrals to calculate areas and volumes. We also give some other typical applications of multivariate
integral calculus.
The relevant MATLAB functions have already been introduced in the previous chapter. They are the following:
syms x, int(f(x), x) or int(‘f(x)’, ‘x’)
Computes the indefinite integral ò f ( x )dx
int (int (‘f(x,y)’, ‘x’), ‘y’)
Computes the double integral òò f ( x , y )dxdy
syms x and int (int (f(x,y), x), y)
Computes the double integral òò f ( x , y )dxdy
int (int (int (… int (‘f(x,y…z)’, ‘x’), ‘y’)…), ‘z’)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz .
syms x y z, int(int(int(…int(f(x,y…z), x), y)…), z)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz
syms x a b, int (f (x), x, a, b)
b
Computes the definite integral ò a
f ( x )dx
int (‘f (x)’, ‘x’, ‘a’, ‘b’)
b
Computes the definite integral ò a
f ( x )dx
int(int(‘f(x,y)’, ‘x’, ‘a’, ‘b’), ‘y’, ‘c’, ‘d’))
b d
Computes the definite integral ò ò
a c
f ( x , y )dxdy
syms x y b c d, int (int (f(x,y), x, a, b), y, c, d)
b d
Computes the definite integral ò ò
a c
f ( x , y )dxdy
163
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
int(int(int(…int(‘f(x,y,…,z)’, ‘x’, ‘a’, ‘b’), ‘y’, ‘c’, ‘d’),…), ‘z’, ‘e’, ‘f’)
b d f
Computes ò ò ¼ò f ( x , y ,¼ , z ) dxdy ¼ dz
a c e
In addition, with prior use of the command maple, the following commands can be used for integration:
int (f (x), x)
Computes the indefinite integral ò f ( x )dx
int (f (x), x = a…b)
b
Computes the definite integral òa
f ( x )dx
int(f(x), x=a..b, continuous)
Computes the definite integral removing the check for continuity of
the function in the interval (a, b)
int (int (f(x,y), x), y)
Computes the indefinite integral òò f ( x , y )dxdy
int(int(int(…int(f(x,y,…,z), x), y),…), z)
Computes ò ò ¼ ò f ( x , y ,¼, z ) dxdy ¼dz
int (int (f(x,y), x = a..b, y = c…d))
b d
Computes the definite integral ò ò
a c
f ( x , y )dxdy
int (int (int (… int (f(x,y,…,z), x = a…b, y = c…d,…, z = e…f))))
b d f
Computes ò ò ¼ò f ( x , y ,¼ , z ) dxdy ¼ dz
a c e
Commands for multivariate integration are also available in the student package, again requiring the prior use of
the command maple. Their syntax is as follows:
Doubleint (expr, var1, var2, name)
Finds the inert double integral of the expression in the variables var1
and var2, where the domain of integration is specified by name
Doubleint (expr, var1=a1..b1, var2=a2..b2)
Finds the specified definite inert double integral
Tripleint (expr, var1, var2, var3, name)
Finds the specified indefinite inert triple integral
Tripleint (expr, var1=a1..b1, var2 = a2…b2, var3 = a3…b3)
Finds the specified inert definite triple integral
164
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
changevar (expr1=expr2,Fnc,var)
Changes the variable in the expression Fnc (Int, Limit, or Sum). The
new variable is defined in terms of the old where expr1 is an expression
in the old variable (appearing in Fnc) and expr2 is an expression in the
new variable var. Used to change variables in integrals.
changevar ({eqn1,eqn2},Doubleint,[var1,var2])
Change of variable in double integrals. The equations defining the new
variables in terms of the old are in the form expr1i = expr2i (i = 1, 2),
where the expr1i are expressions in the old variables, and expr2i are
expressions in the new variables var1 and var2.
changevar ({eqn1,eqn2,eqn3},Tripleint,[var1,var2,var3])
Change of variables for triple integrals
completesquare (expression)
Completes the square with respect to all variables appearing in the
specified polynomial expression
completesquare (expr, variable)
Completes the square of the specified polynomial expression with
respect to the given variable
integrand (expression)
Extracts all integrands of an expression containing integrals (Int,
Doubleint, Tripleint)
powsubs (expr1=expr2,expr)
Replaces all occurrences of expr1 in expr by expr2 (expr1 is a
subexpression of expr)
powsubs ({eqn1,…,eqnn},expr) or subs([eqn1,…,eqnn],expr)
Simultaneously performs the substitutions in expr described by the n
equations eqni of the form expr1i = expr2i, i = 1… n
A = òò dx dy
S
If, for example, S is determined by a < x < b and f (x) < y < g(x), the area will be:
b g (x)
A = ò dx ò dy
a f (x)
165
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
d h(b)
A = ò dy ò dx
c h(a )
If the region S is determined by curves whose equations are given in polar coordinates with radius vector r and
angle a, its area A is given by the formula:
A = òò r da dr
S
t g (a )
A = ò da ò r dr
s f (a )
Exercise 7-1
Calculate the area of the region above the OX axis bounded by the OX axis, the parabola y 2 = 4x and the line
x + y = 3.
First, we create a graphical representation of the problem, which is presented in Figure 7-1:
>> fplot ('[(4*x) ^(1/2), 3-x]', [0,4,0,4])
Figure 7-1.
166
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
We see that in the enclosed region y is limited between 0 and 2 (0 < y < 2) and x is limited between the curves
x = ( y ^ 2)/4 and x = 3 − y. We calculate the requested area as follows:
>> int(int('1','x','y^2/4','3-y'),'y',0,2)
ans =
10/3
Exercise 7-2
Calculate the area in the first quadrant bounded between the semicubical parabola y 2 = x 3 and the bisector of the
first quadrant.
First, we graphically represent the problem (see Figure 7-2):
>> fplot ('[^(3/2) x, x]', [-4.5-1-4])
Figure 7-2.
167
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Inspecting the graph we see that the enclosed area is defined for x ranging between 0 and 1 (0 < x < 1) and y
ranging between the curves y = x and y = x ^( 3/2). Therefore the area of the requested region can be calculated as:
>> int(int('1','y','x^(3/2)','x'),'x',0,1)
ans =
1/10
Exercise 7-3
Calculate the area outside the circle with polar equation r = 2, and inside the cardioid with polar equation
r = 2 (1 + cos (a )).
First, we graphically represent the problem (see Figure 7-3):
>> a = 0:0. 1:2 * pi;
>> r = 2 * (1 + cos (a));
>> polar (a, r)
>> hold on;
>> r=2*ones(size(a));
>> polar (a, r)
Figure 7-3.
168
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Inspecting the graph, we see that, by symmetry, we can calculate half of the required area by varying a between
0 and Pi / 2 (0 < a < Pi /2) and r between the curves r = 2 and r = 2 (1 + cos (a)):
>> pretty(int(int('r','r',2,'2*(1+cos(a))'),'a',0,pi/2))
1/2 pi + 4
The required area is therefore p + 8 square units.
2 2
æ¶z ö æ¶z ö
S = òò 1+ ç ÷ +ç ÷ dx dy
R
è¶x ø è¶ y ø
The surface area of a surface S defined by x = f(y, z), where (y, z) ranges over a region R in the OYZ plane,
is given by:
2 2
æ¶x ö æ¶x ö
S = òò 1+ ç ÷ +ç ÷ dy dz
R
è¶y ø è¶z ø
The surface area of a surface S defined by y = f (x, z), where (x, z) ranges over a region R in the OXZ plane, is
given by:
2 2
æ¶ y ö æ¶ y ö
S = òò 1+ ç ÷ +ç ÷ dx dz
R
è¶x ø è¶z ø
Exercise 7-4
Calculate the area of the surface of the cone x 2 + y 2 = z 2, limited above the OXY plane and cut by the cylinder
x 2 + y 2 = b y.
The projection of the surface onto the OXY plane is the disc bounded by the circle with equation x 2+ y 2 = by.
Then we can find the surface area via the first formula above as follows:
>> maple('z:=(x,y)->sqrt((x^2+y^2)/a)');
>> pretty(int(int('(1+diff(z(x,y),x)^2+diff(z(x,y),y)^2)^(1/2)','x','-(b*y-y^2)^(1/2)',
'(b*y-y^2)^(1/2)'),'y',0,'b'))
2 1/2
b pi (a + 1)
1/4 -------------
1/2
a
169
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-5
Calculate the area of the paraboloid x 2 + y 2 = 2 z limited below the plane z = 2.
>> maple('z:=(x,y)->(x^2+y^2)/2');
>> pretty(maple('(1+diff(z(x,y),x)^2+diff(z(x,y),y)^2)^(1/2)'))
2 2 1/2
(1 + x + y )
This expression is suitable for a transformation to polar coordinates:
>> maple('m:=(x,y)->sqrt(1+x^2+y^2)');
>> pretty(simple('m(r*cos(a),r*sin(a))'))
2 1/2
(1 + r )
The requested integral is calculated as 4 times the integral delimited by the first quadrant of the circle r = 2 :
>> pretty(simple(int(int('r*sqrt(1+r^2)','r',0,2),'a',0,pi/2)))
1/2
1/6 (5 5 - 1) pi
( )
The result of the integral will be 4 times the previous value, i.e. 2 5 5 - 1 p/ 3 .
Figure 7-4 shows the graphical representation of the problem:
>> [x, y] = meshgrid(-3:.1:3);
>> z=(1/2)*(x.^2+y.^2);
>> mesh(x,y,z)
>> hold on;
>> z=2*ones(size(z));
>> mesh(x,y,z)
>> view(-10,10)
170
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Figure 7-4.
V = òò f ( x , y ) dx dy = òò z dx dy
R R
The volume V of a cylindroid limited in its upper part by the surface with equation x = f(y, z), its lower part by the
OYZ plane and laterally by the straight cylindrical surface which cuts the OYZ plane bordering a region R, is:
V = òò f ( y , z ) dy dz = òò x dy dz
R R
The volume V of a cylindroid limited in its upper part by the surface with equation y = f(x, z), its lower part by the
OXZ plane and laterally by the straight cylindrical surface which cuts the OXZ plane bordering a region R, is:
V = òò f ( x , z ) dx dz = òò y dx dz
R R
171
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-6
Calculate the volume in the first octant bounded between the OXY plane, the plane z = x + y + 2 and the cylinder
x 2 + y 2 = 16.
We begin by graphically representing the problem (see Figure 7-5) with the plane in cartesian coordinates and
parameterizing the cylinder:
>> t =(0:.1:2*pi);
>> u =(0:.1:10);
>> x = 4 * cos (t)'* ones (size (u));
>> y = 4 * sin (t)'* ones (size (u));
>> z = ones (size (t))'* u;
>> mesh (x, y, z)
>> hold on;
>> [x,y]=meshgrid(-4:.1:4);
>> z = x + y + 2;
>> mesh (x, y, z)
>> set(gca,'Box','on');
>> view(15,45)
Figure 7-5.
172
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-7
Calculate the volume bounded by the paraboloid x 2 + 4 y 2 = z and laterally by the cylinders with equations y 2 = x
and x 2 = y.
First, we graphically represent the problem (see Figure 7-6):
>> [x,y]=meshgrid(-1/2:.02:1/2,-1/4:.01:1/4);
>> z = x ^ 2 + 4 * y. ^ 2;
>> mesh(x,y,z)
>> hold on;
>> y=x.^2;
>> mesh(x,y,z)
>> hold on;
>> x=y.^2;
>> mesh(x,y,z)
>> set(gca,'Box','on')
>> view(-60,40)
Figure 7-6.
173
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
òòò dx dy dz
R
The volume of a three-dimensional region R whose boundary surface equations are expressed in cylindrical
coordinates is given by the triple integral:
òòò r dz dr da
R
The volume of a three-dimensional region R whose boundary surfaces equations are expressed in spherical
coordinates is given by the triple integral:
òòò r
2
sin ( b ) drdbda
R
Exercise 7-8
Calculate the volume bounded by the paraboloid x 2 + y 2 = z and the cylinder with equation z = a 2 − y 2.
The volume will be four times the following integral:
>> pretty(simple(int(int(int('1','z','a*x^2+y^2','a^2-y^2'),'y',0,
'sqrt((a^2-a*x^2)/2)'),'x',0,'sqrt(a)')))
4
a (log(-a) - log(a))
1/8 ---------------------
1/2
(-2 a)
174
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
We can assume a > 0 and simplify the value of 4 times the previous integral as follows (note this gives the real
part of the given expression):
>> pretty(simple(4*int(int(int('1','z','a*x^2+y^2','a^2-y^2'),'y',0,
'sqrt((a^2-a*x^2)/2)'),'x',0,'sqrt(a)')))
7/2 1/2
1/4 a 2 pi
Exercise 7-9
Figure 7-7.
175
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-10
Calculate the volume enclosed by the cylinder r = 4cos(a), the plane z = 0 and the sphere with equation r 2 + z 2 = 16.
The volume of the enclosure, given in cylindrical coordinates, is:
>> pretty(simple(int(int(int('r','z',0,'sqrt(16-r^2)'),'r',0,'4*cos(a)'),'a',0,pi)))
256/9 + 64/3 pi
Exercise 7-11
Calculate the volume enclosed by the cone b = p /4 and the sphere with equation r = 2 k cos (b).
The volume of the enclosure, given in spherical coordinates, will be four times the result of the following integral:
>> pretty(simple(int( int(int('r^2*sin(b)','r',0,'2*k*cos(b)'),'b',0,pi/4),'a',0,pi/2)))
3
1/4 k pi
The final result is therefore k3p.
æ¶n ¶m ö
ò m( x , y ) dx + n( x , y ) dy = òò çè ¶ x - ¶ y ÷ø dA
C R
176
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-12
ò (x +e ) dx + (2 y + cos(x ))dy
y
C
where C is the boundary of the region enclosed by the parabolas y = x 2 , x = y 2.
The two parabolas intersect at the points (0,0) and (1,1). We graphically represent the problem (see Figure 7-8)
and calculate the integral:
plot([x^2,sqrt(x)],x,0,1.2);
Figure 7-8.
-.67644120041679251512532326651204
177
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
òòS
f • n dS = òòò Div( f ) dV
Q
The left-hand side of this equality is called the outflow of the vector field f across the surface S.
Exercise 7-13
Use the divergence theorem to find the outflow of the vector field f = (xy + x 2, yz + xy 2, z z, xz + xyz 2 ) across the
surface of the cube in the first octant bounded by the planes x = 2, y = 2 and z = 2.
>> maple ('f: = vector([x*y+x^2*z*y, y*z+x*y^2*z, x*z+x*y*z^2])');
>> maple ('v: = vector ([x, y, z])');
>> int(int(int('diverge(f,v)','x',0,2),'y',0,2),'z',0,2)
ans =
72
ò F × dr = ò ò ( curl F ) × n ds
C S
178
www.pdfgrip.com
Chapter 7 ■ Integration in Several Variables and Applications
Exercise 7-14
òC
- y 3dx + x 3dy - z 3dz
where C is the intersection of the cylinder x 2 + y 2 = 1 and the plane x + y + z = 1, and the orientation of C
corresponds to counterclockwise rotation of the OXY plane.
The curve C bounds the surface S defined by z = 1 − x − y = f (x,y) for (x, y) in the domain D = {(x,y) / x 2 + y 2 = 1}.
We put F = − y 3 i + x 3 j − z 3 k.
Now, we calculate the curl of F and integrate over the surface S:
>> maple ('F: = vector([-y^3,x^3,z^3])');
>> maple ('v: = vector ([x, y, z])');
>> pretty(sym(maple('curl(F,v)')))
2 2
[0, 0, 3 x + 3 y]
Therefore, we have to calculate the integral ∫D (3x 2 + 3y 2 ) dx dy. Changing to polar coordinates, we obtain:
>> pretty(simple(int(int('3*r^3', 'a',0,2*pi),'r',0,1)))
3/2 pi
179
www.pdfgrip.com
Chapter 8
Differential Equations
Although it implements only a relatively small number of commands related to this topic, MATLAB’s treatment of
differential equations is nevertheless very efficient. We shall see how we can use these commands to solve each type
of differential equation algebraically. Numerical methods for the approximate solution of equations and systems of
equations are also implemented.
The basic command used to solve differential equations is dsolve. This command finds symbolic solutions
of ordinary differential equations and systems of ordinary differential equations. The equations are specified by
symbolic expressions where the letter D is used to denote differentiation, or D2, D3, etc, to denote differentiation of
order 2,3,..., etc. The letter preceded by D (or D2, etc) is the dependent variable (which is usually y), and any letter
that is not preceded by D (or D2, etc) is a candidate for the independent variable. If the independent variable is not
specified, it is taken to be x by default. If x is specified as the dependent variable, then the independent variable is t.
That is, x is the independent variable by default, unless it is declared as the dependent variable, in which case the
independent variable is understood to be t.
You can specify initial conditions using additional equations, which take the form y(a) = b or Dy(a) = b,..., etc.
If the initial conditions are not specified, the solutions of the differential equations will contain constants of integration,
C1, C2,..., etc. The most important MATLAB commands that solve differential equations are the following:
dsolve( ‘equation’, ‘v’)
Solves the given differential equation, where v is the independent variable (if ‘v’ is
not specified, the independent variable is x by default). This returns only explicit
solutions.
dsolve( ‘equation’, ‘condition_initial’,..., ‘v’)
Solves the given differential equation subject to the specified initial condition
dsolve( ‘equation’, ‘cond1’, ‘cond2’,..., ‘condn’, ‘v’)
Solves the given differential equation subject to the specified initial conditions
dsolve( ‘equation’, ‘cond1,cond2,...,condn’, ‘v’)
Solves the given differential equation subject to the specified initial conditions
dsolve( ‘eq1’, ‘eq2’,..., ‘eqn’, ‘cond1’, ‘cond2’,..., ‘condn’, ‘v’)
Solves the given differential system subject to the specified initial conditions
(explicit)
dsolve( ‘eq1,eq2,...,eqn’, ‘cond1,cond2,...,condn’ , ‘v’)
Solves the given differential system subject to the specified initial conditions
181
www.pdfgrip.com
Chapter 8 ■ Differential Equations
There is another group of important MATLAB commands that solve differential equations, requiring the prior use
of the command maple. They are as follows:
dsolve (eqn, fnc (var))
Symbolically solves the ordinary differential equation eqn for the function
fnc(var). The solution is usually returned in implicit form as an equation in var,
fnc (var) and the implied constants C_1,..., C_n. var is an independent variable
and fnc is the dependent variable. The result may also contain calls to the DEsol
command.
dsolve (expr, fnc (var))
Solves the differential equation expr = 0
dsolve({deqn,cond1,...,condn},fnc(var))
Solves the differential equation deqn subject to the initial conditions
cond1,..., condn
dsolve (deqn, fnc (var), explicit = true)
Gives the solution in explicit form of the differential equation deqn, if possible
dsolve(deqn,fnc(var),method=laplace)
Forces the use of the Laplace transform method in the solution of the differential
equation
dsolve(deqn,fnc(var),type=series)
Forces the use of the Taylor series method in the solution of the differential
equation, with the degree of the polynomial solution determined by the global
variable Order
dsolve(deqn,fnc(var),output=basis)
The solution is given in terms of the specified basis
powseries[powsolve](deqn,cond1,...,condn)
Solves the differential equation deqn subject to the initial conditions given in
terms of power series
dsolve ({deqn1,..., deqnn}, {(var) fnc1... fncn (var)})
Solves the specified system of differential equations for the independent variable
var and the dependent variables fnc1,..., fncn
dsolve ({deqn1,..., deqnn, cond1,..., condm}, {(var) fnc1..)fncn (var)})
Solves the specified system of differential equations for the independent variable
var and the dependent variables fnc1,..., fncn, subject to the initial conditions
cond1,..., condm
DEsol (deqn, fnc (var))
Gives the solution of the differential equation deqn for the function fnc(var) in
inert form
182
www.pdfgrip.com
Chapter 8 ■ Differential Equations
183
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Next we solve a couple of systems, both with and without initial values:
>> pretty(dsolve('Dx = y', 'Dy = -x'))
x (t) = (t) C1 + C2 cos (t), y (t) = C1 cos (t) - C2 sin (t)
>> S=dsolve('Df = 3*f+4*g', 'Dg = -4*f+3*g');[S.f, S.g]
ans =
[exp(3*t)*cos(4*t)*C1+exp(3*t)*sin(4*t)*C2, -exp(3*t)*sin(4*t)*C1+ exp(3*t)*cos(4*t)*C2]
>> T= dsolve('Df = 3*f+4*g, Dg = -4*f+3*g', 'f(0)=0, g(0)=1') ;[T.f, T.g]
ans =
[ exp(3*t)*sin(4*t), exp(3*t)*cos(4*t)]
This last system can also be solved in the following way:
>> pretty(maple('dsolve({diff(f(x),x)= 3*f(x)+4*g(x), diff(g(x),x)=- 4*f(x) + 3*g(x),
f(0)=0,g(0)=1}, {f(x),g(x)})'))
{f (x) = exp (3 x) sin(4 x), g (x) = exp(3 x) cos(4 x)}
We now solve a differential equation by the Laplace transform method:
>> dsolve ({diff (y (t), t$ 2) + 5 * diff (y (t), t) + 6 * y (t) = 0, y (y) (0) (0) = 0, D = 1},
y (t), method = laplace);
y (t) = - exp(-3 t) + exp(-2 t)
Next we solve a differential system in series form:
>> dsolve({diff(y(x),x)=z(x)-y(x)-x,diff(z(x),x)=y(x),y(0)=0,z(0)=1},{y(x), z(x)}, type=series);
2 3 4 5 6
{y(x) = x - x + 1/2 x - 5/24 x + 1/15 x + O(x ),
2 3 4 5 6
z(x) = 1 + 1/2 x - 1/3 x + 1/8 x - 1/24 x + O(x )}
184
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-1
185
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-2
186
www.pdfgrip.com
Chapter 8 ■ Differential Equations
>> pretty(sym(maple('collect(v^2*y^3*d(v)+v^3*y^2*d(y)-y^3*d(v),{d(v),d(y)})')))
3 2 2 3 3
v y d(y) + (v y - y ) d(v)
If we divide the previous expression by v 3y 3, and group the terms in d (v) and d (y), we already have an equation in
separated variables.
>> pretty(sym(maple('collect(((v^2*y^3-y^3)*d(v)+v^3*y^2*d(y))/(v^3*y^3), {d(v),d(y)})')))
2 3 3
(v y - y ) d(v) d(y)
----------------- + ----
3 3 y
v y
The previous expression can be simplified.
>> pretty(maple('convert(collect(((v^2*y^3-y^3)*d(v)+v^3*y^2*d(y))/(v^3*y^3),
{d(v),d(y)}),parfrac,y)'))
2
(v - 1) d(v) d(y)
------------- + ----
3 y
v
Now, we solve the equation:
>> pretty(simple(sym('int((v^2-1)/v^3,v)+int(1/y,y)')))
1
log(v) + ---- + log(y)
2
2 v
Finally we reverse the change of variable:
>> pretty(simple(sym('subs(v=x/y,log(v)+1/2/v^2+log(y))')))
2
y
log(x) + 1/2 ----
2
x
Thus the general solution of the original differential equation is:
1 y2
log ( x ) + = C.
2 x2
187
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Now we can represent the solutions of this differential equation graphically. To do this we graph the solutions with
parameter C, which is equivalent to the following contour plot of the function defined by the left-hand side of the
above general solution (see Figure 8-1):
>> [x,y]=meshgrid(0.1:0.05:1/2,-1:0.05:1);
>> z=y.^2./(2*x.^2)+log(x);
>> contour(z,65)
Figure 8-1.
is said to be exact if ∂N/∂x = ∂M/∂y. If the equation is exact, then there exists a function F such that its total differential
dF coincides with the left-hand side of the above equation, i.e.:
dF = M(x,y) dx + N(x,y) dy
therefore the family of solutions is given by F(x,y) = C.
The exercise below follows the usual steps of an algebraic solution to this type of equation.
188
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-3
189
www.pdfgrip.com
Chapter 8 ■ Differential Equations
>> pretty (simple (sym ('solve (x * exp(y*x) + x * cos(y*x) + diff (g (y), y) = n(x,y),))))
((((diff (g (y), y))')))
ans = 1
Thus g'(y) = 1, so the final solution will be, omitting the addition of a constant:
>> pretty (simple (sym (maple ('subs (g (y) = int(1,y),-x+exp(y*x) + sin(y*x) + g (y))'))))
-x + exp(y x) + sin(y x) + y
To graphically represent the family of solutions, we draw the following contour plot of the above expression
(Figure 8-2):
>> [x,y]=meshgrid(-2*pi/3:.2:2*pi/3);
>> z =-x+exp (y.*x) + free (y.*x) + y;
>> contour(z,100)
Figure 8-2.
In the following section we will see how any reducible differential equation can be transformed to an exact
equation using an integrating factor.
190
www.pdfgrip.com
Chapter 8 ■ Differential Equations
eò
P ( x ) dx
MATLAB implements these solutions of linear differential equations, and offers them whenever the integral
appearing in the integrating factor can be found.
Exercise 8-4
= f ( x ).
If the function f(x) is identically zero, the equation is called homogeneous. Otherwise, the equation is called
non-homogeneous. If the functions ai(x) (i = 1, ..., n) are constant, the equation is said to have constant coefficients.
191
www.pdfgrip.com
Chapter 8 ■ Differential Equations
A concept of great importance in this context is that of a set of linearly independent functions. A set of functions
{f1(x), f2(x), ..., fn(x)} is linearly independent if, for any x in their common domain of definition, the Wronskian
determinant of the functions is non-zero. The Wronskian determinant of the given set of functions, at a point x of
their common domain of definition, is defined as follows:
f1 ( x ) f2 (x ) f3 (x ) ¼ fn ( x )
f1¢ ( x ) f 2¢ ( x ) f 3¢ ( x ) ¼ f n¢ ( x )
f1¢¢ ( x ) f 2¢¢ ( x ) f3 (x ) ¼
¢¢
f n¢¢ ( x ) = W (x)
¼ ¼ ¼ ¼ ¼
f1(n-1) ( x ) f 2(n-1) ( x ) f 3(n-1) ( x ) ¼ f n(n-1) ( x )
The MATLAB command maple('Wronskian') allows you to calculate the Wronskian matrix of a set of functions.
Its syntax is:
maple ('Wronskian(V,variable)')
computes the Wronskian matrix corresponding to the vector of functions V with
independent variable x.
A set S = {f1(x), ..., fn(x)} of linearly independent non-trivial solutions of a homogeneous linear equation of order n
a0 ( x ) y ( x ) + a1 ( x ) y ¢ ( x ) + a2 ( x ) y ¢¢ ( x ) +¼+ an ( x ) y (n ) ( x ) = 0
is called the characteristic equation of the homogeneous differential equation with constant coefficients.
The solutions of this characteristic equation determine the general solutions of the corresponding differential
equation.
192
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-5
is linearly independent.
>> W=sym(maple('Wronskian(vector([exp(x),x*exp(x), x^2*exp(x)]),x)'))
W =
[exp(x), x*exp(x), x^2*exp(x)]
[exp(x), exp(x)+x*exp(x), 2*x*exp(x)+x^2*exp(x)]
[exp(x), 2*exp(x)+x*exp(x), 2*exp(x)+4*x*exp(x)+x^2*exp(x)]
>> pretty(determ(W))
3
2 exp(x)
This gives us the value of the Wronskian, which is obviously always non-zero. Therefore the set of functions is
linearly independent.
åa ( x ) y ( ) ( x ) = a ( x ) y ( x ) + a ( x ) y¢ ( x ) + a ( x ) y ² ( x ) + ¼ + a ( x ) y ( ) ( x )
k =0
k
k
0 1 2 n
n
=0
is said to have constant coefficients if the functions ai(x) (i = 1, ..., n) are all constant (i.e. they do not depend on the
variable x).
The equation:
n
a0 + a1m + a2m 2 + ¼ + anm n = åaim i = 0
i =0
is called the characteristic equation of the above differential equation. The solutions (m1, m2, ..., mn) of this
characteristic equation determine the general solution of the associated differential equation.
If the mi (i = 1, ..., n) are all different, the general solution of the homogeneous equation with constant coefficients is:
y ( x ) = c1e m1 x + c 2e m2 x + ¼ + cn e mn x
ci e mi x + ci +1 xe mi x + ci +2 x 2e mi x + ¼ + ci +k x k e mi x .
193
www.pdfgrip.com
Chapter 8 ■ Differential Equations
If the characteristic equation has a complex root mj = a + bi , then its complex conjugate mj +1 = a − bi is also a
root. These two roots determine a pair of terms in the general solution of the homogeneous equation:
MATLAB directly applies this method to obtain the solutions of homogeneous linear equations with constant
coefficients, using the command dsolve or maple('dsolve').
Exercise 8-6
2 y ² + 5 y ¢ + 5 y = 0, y (0) = 0, y ¢(0) = 1 2
>> pretty(dsolve('3*D2y+2*Dy-5*y=0'))
C1 exp(t) + C2 exp(- 5/3 t)
>> pretty(dsolve('2*D2y+5*Dy+5*y=0','y(0)=0,Dy(0)=1/2'))
1/2 1/2 1/2 1/2
2/15 3 5 exp(- 5/4 t) sin(1/4 3 5 t)
Exercise 8-7
194
www.pdfgrip.com
Chapter 8 ■ Differential Equations
åa ( x ) y ( ) ( x ) = a ( x ) y ( x ) + a ( x ) y¢ ( x ) + a ( x ) y ² ( x ) + ¼ + a ( x ) y ( ) ( x )
k =0
k
k
0 1 2 n
n
= f ( x ).
Suppose {y1(x), y2(x),....., yn(x)} is a linearly independent set of solutions of the corresponding homogeneous
equation:
a0 ( x ) y ( x ) + a1 ( x ) y ¢ ( x ) + a2 ( x ) y ² ( x ) + ¼ + an ( x ) y (n ) ( x ) = 0.
f ( x )Wi [ y1 ( x ) , y 2 ( x ) ,¼, yn ( x )]
ui ( x ) = ò dx ( i = 1,¼,n ) .
W [ y1 ( x ) , y 2 ( x ) ,¼, yn ( x )]
Here Wi[y1(x), y2(x), ..., yn(x)] is the determinant of the matrix obtained by replacing the i-th column of the
Wronskian matrix W[y1(x), y2(x), ..., yn(x)] by the transpose of the vector (0,0,..., 0, 1).
The solution of the non-homogeneous equation is then given by combining the general solution of the
homogeneous equation with the particular solution of the non-homogeneous equation. If the roots mi of the
characteristic equation of the homogeneous equation are all different, the general solution of the non-homogeneous
equation is:
y ( x ) = c1e m1 x + c 2e m2 x + ¼ + cn e mn x + y p ( x ) .
If some of the roots are repeated, we refer to the general form of the solution of a homogeneous equation
discussed earlier.
195
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-8
We will follow the algebraic method of variation of parameters to solve the first equation. We first consider the
characteristic equation of the homogeneous equation to obtain a set of linearly independent solutions.
>> solve('m^2+4*m+13=0')
ans =
[- 2 + 3 * i]
[- 2 - 3 * i]
>> maple ('f: = x - > x * cos(3*x) ^ 2');
>> maple ('y1: = x - > exp(-2*x) * cos(3*x)');
>> maple ('y2: = x - > exp(-2*x) * sin(3*x)');
>> maple ('W: = x - > Wronskian ([y1 (x), y2 (x)], x)');
>> pretty(simple(sym(maple('det(W(x))'))))
3 exp(-4 x)
We see that the Wronskian is non-zero, indicating that the functions are linearly independent. Now we calculate
the functions Wi (x ), i = 1, 2.
>> maple ('W1: x-= > array ([[0, y2 (x)], [1, diff ((y2) (x), x)]])');
>> pretty(simple(sym(maple('det(W1(x))'))))
-exp(-2 x) sin (3 x)
>> maple ('W2: x-= > array ([[y1 (x), 0], [diff ((y1) (x), x), 1]])');
>> pretty(simple(sym(maple('det(W2(x))'))))
exp(-2 x) cos(3 x)
Now we calculate the particular solution of the non-homogeneous equation.
>> maple('u1:=x->factor(simplify(int(f(x)*det(W1(x))/det(W(x)),x)))');
>> maple ('u1 (x)')
ans =
1/14652300*exp(2*x)*(129285*cos(9*x)*x-6084*cos(9*x)-28730*sin(9*x)*x-
13013*sin(9*x)+281775*cos(3*x)*x-86700*cos(3*x)-187850*sin(3*x)*x-36125*sin(3*x))
196
www.pdfgrip.com
Chapter 8 ■ Differential Equations
>> maple('u2:=x->factor(simplify(int(f(x)*det(W2(x))/det(W(x)),x)))');
>> maple ('u2 (x)')
ans =
1/14652300 * exp(2*x) * (563550 * cos(3*x) * x+108375 * cos(3*x) + 845325 * sin(3*x) *
x-260100 * sin(3*x) + 28730 * cos(9*x) * x+13013 * cos(9*x) + 129285 * sin(9*x) * x-6084 *
sin(9*x))
>> maple ('yp: = x - > factor (simplify (y1 (x) * (x) u1 + y2 (x) * u2 (x)))');
>> maple('yp(x)')
ans =
-23/1105 * x * cos(3*x) ^ 2 + 13436/1221025 * cos(3*x) ^ 2 + 24/1105 * cos(3*x) * sin(3*x) *
x + 3852/1221025 * cos(3*x) * sin(3*x) + 54/1105 * x-21168/1221025
Then we can write the general solution of the non-homogeneous equation:
>> maple (' y: = x - > simplify (c1 * y1 (x) + c2 * y2 (x) + yp (x))');
>> maple ('combine (y (x), trig)')
ans =
C1 * exp(-2*x) * cos(3*x) + c2 * exp(-2*x) * sin(3*x)-23/2210 * x * cos(6*x) + 1/26 * x +
6718/1221025 * cos(6*x)-2/169 + 12/1105 * x * sin(6*x) + 1926/1221025 * sin(6*x)
Now we graphically represent a set of solutions, for certain values of c1 and c2 (see Figure 8-3)
>> ezplot(simple(sym(maple('subs(c1=-5,c2=4,y(x))'))),[-1,1])
>> hold on
>> ezplot(simple(sym(maple('subs(c1=-5,c2=-4,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=-5,c2=-2,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=-5,c2=2,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=-5,c2=1,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=-5,c2=-1,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=-1,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=1,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=2,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=-2,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=4,y(x))'))),[-1,1])
>> ezplot(simple(sym(maple('subs(c1=5,c2=-4,y(x))'))),[-1,1])
197
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Figure 8-3.
There is a MATLAB command that allows you to solve differential equations by the method of variation of
parameters directly. Its syntax is as follows:
maple ('DEtools [varparam]([expr1,...,exprn], exprrhs, varind)')
Finds the general solution by the method of variation of parameters of an ordinary
differential equation in the independent variable varind, where exprrhs is the right-hand
side of the original equation (non-homogeneous part) and the expri (i = 1... n) are the
solutions of the corresponding homogeneous equation of order n. The command value can
help to obtain the final solution.
For the second differential equation we directly apply dsolve, obtaining the solution. First, we find the general
solution of the homogeneous equation:
>> maple ('solhomog: = dsolve (diff (y (x), x$ 2) - 2 * diff (y (x), x) + y (x) = 0 y (x))')
ans =
solhomog: = y(x) = _C1 * exp (x) + _C2 * exp (x) * x
As the differential equation is of order 2 we must take two solutions (for two pairs of any values of C1 and C2).
The simplest solutions that can be taken are those which form a basis for the set of solutions (exp (x) and x exp (x)),
which can be found directly as follows:
>> maple('solhomog:=dsolve(diff(y(x),x$2)-2*diff(y(x),x)+y(x)=0,y(x),output=basis)')
ans =
solhomog := [exp(x), exp(x)*x]
198
www.pdfgrip.com
Chapter 8 ■ Differential Equations
As the inhomogeneous part of the right-hand side of the equation is e x ln (x ), and the independent variable is x,
we can already apply the command varparam as follows:
>> pretty(sym(maple('DEtools[varparam](solhomog,exp(x)*ln(x),x)')))
2 2
_C[1] exp(x) + _C[2] exp(x) x + (- 1/2 x log(x) + 1/4 x ) exp(x)
+ (x log(x) - x) exp(x) x
The second differential equation is solved directly as follows:
>> pretty(simple(dsolve('D2y-2*Dy+y=exp(t)*log(t)')))
2 2
1/4 exp(t) (2 t log(t) - 3 t + 4 C1 + 4 C2 t)
åa x
k =0
k
k
y ( k ) ( x ) = a0 y ( x ) + a1 xy ¢ ( x ) + a2 x 2 y ² ( x ) +¼+ an x n y (n ) ( x )
= f (x)
is called a Cauchy–Euler equation.
This equation can be reduced to a homogeneous linear equation with constant coefficients by replacing x = et.
This leads us to solve the equation:
MATLAB solves this type of equation directly with the command dsolve or maple('dsolve').
199
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-9
MATLAB provides the commands maple('laplace') and maple('invlaplace') to calculate the Laplace transform
and inverse Laplace transform of an expression with respect to a variable. The syntax is as follows:
maple ('laplace (expression, t, s)')
Calculates the Laplace transform of a given expression with respect to t.
The transformed variable is s.
maple ('invlaplace (expression, s, t)')
Computes the inverse Laplace transform of the given expression with respect to s.
The inverse variable is t.
Here are some examples:
>> pretty(sym(maple('laplace(t^(3/2)-exp(t)+sinh(a*t), t, s)')))
1/2
pi 1 a
3/4 ----- - ----- + -------
5/2 s - 1 2 2
s s - a
>> pretty(sym(maple('invlaplace(s^2/(s^2+a^2)^(3/2), s, t)')))
- t BesselJ(1, a t) a + BesselJ(0, a t)
The Laplace transform and its inverse are used to solve certain differential equations. The method is to calculate
the Laplace transform of each term of the equation to obtain a new differential equation, which we then solve. Finally,
we find the solution of the original equation by applying the inverse Laplace transform to the solutions just found.
200
www.pdfgrip.com
Chapter 8 ■ Differential Equations
MATLAB provides the ‘laplace’ option in the maple(‘dsolve’) command, which forces the program to solve the
differential equation using the Laplace transform method. The syntax is as follows:
maple ('dsolve (equation, func (var), 'laplace'))
Exercise 8-10
201
www.pdfgrip.com
Chapter 8 ■ Differential Equations
This gives the solution of the Laplace transformed equation. To calculate the solution of the original equation we
calculate the inverse Laplace transform of the solution obtained in the previous step.
>> maple('TL0:=s->simplify(subs(y(0)=1,(D(y))(0)=1,TL(s)))');
>> solution=simple(sym(maple('invlaplace(TL0(s),s,x)')));
>> pretty(solution)
1/2 1/2 1/2
1 sin(3 x) 3 35 cos(3 x)
1/4 x - 1/8 - -------- + 5/8 ---------------- + ---- -----------
3 exp(x) exp(x) 24 exp(x)
This gives the solution of the original differential equation.
We could also have solved it directly via:
>> pretty(simple(sym(maple('dsolve({diff(y(x),x$2)+2*diff(y(x),x)+4*y(x)= x-exp(-x),
y(0)=1,D(y)(0)=1},y(x),laplace)'))))
1/2 1/2 1/2
1 sin(3 x) 3 35 cos(3 x)
y(x) = 1/4 x - 1/8 - -------- + 5/8 ---------------- + ---- ------------
3 exp(x) exp(x) 24 exp(x)
where the eigenvalues {lk} (k = 1, 2, ... n) corresponding to the eigenvectors {Vk} of the matrix A of the system are all
assumed to be different.
If an eigenvalue lk is a complex number ak + bki, then it generates the following component of the overall solution:
C k1Wk1 e lk t + C k2Wk2 e lk t
where:
Wk1 =
1
2
( ) i
(
Vk + V k cos ( bk t ) + Vk - V k sin ( bkt ) ,
2
)
Wk2 =
i
2
( ) 1
(
Vk - V k cos ( bkt ) - Vk + V k sin ( bkt ) .
2
)
Here Vk is the eigenvector corresponding to the eigenvalue lk and V k is its conjugate.
If there is an eigenvalue li of multiplicity m > 1, then it will generate a portion of the general solution of the form:
MATLAB can solve this type of system directly, simply by using the command dsolve or maple('dsolve') with the
familiar syntax.
202
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-11
y ¢ = -2 x + 10 y .
>> S = dsolve('Dx=-5*x+3*y,Dy=-2*x-10*y','t');
>> pretty(sym([S.x,S.y]))
[-2 C1 exp(-8 t) + 3 C1 exp(-7 t) + 3 C2 exp(-7 t) - 3 C2 exp(-8 t) ,
-2 C1 exp(-7 t) + 2 C1 exp(-8 t) + 3 C2 exp(-8 t) - 2 C2 exp(-7 t)]
You can also use the following syntax:
>> pretty (sym (maple ('dsolve ({diff (x (t), t) =-5 * x (t) + 3 * y (t), diff (y (t), t) =
-2 * x (t) - 10 * y (t)}, {x (t), y (t)})')))
{x(t) = -2 _C1 exp(-8 t) + 3 _C1 exp(-7 t) + 3 _C2 exp(-7 t) - 3 _C2 exp(-8 t)
y(t) = -2 _C1 exp(-7 t) + 2 _C1 exp(-8 t) + 3 _C2 exp(-8 t) - 2 _C2 exp(-7 t)}
The general solution of the non-homogeneous system will be X = F(t )C + Xp, which is, using the previous
expression:
X = F ( t )C + F(t )òF -1 ( t ) F ( t ) dt .
This method is a generalization to systems of equations of the method of variation of parameters for simple
equations.
MATLAB can solve such systems of equations directly with the command dsolve or maple('dsolve'), provided the
integrals that appear in the solution can be evaluated.
203
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-12
204
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-13
4 x 2 y ¢¢ + 4 xy ¢ + ( x 2 - 1) y = 0 ,
>> pretty(simple(sym(maple('convert(",polynom)'))))
2 4
1/2 2 4
_C2 (1 - 1/8 x + 1/384 x )
y(x) = _C1 x (1 - 1/24 x + 1/1920 x ) + ---------------------------
1/2
x
>> pretty(simple(sym(maple('dsolve({y(x)*diff(y(x),x$2)+diff(y(x),x)^2+1=0,
y(0)=1,D(y)(0)=1},y(x),series)'))))
2 3 4 5 6
y(x) = 1 + x - x + x - 3/2 x + 5/2 x + O(x )
Exercise 8-14
Solve the following two systems of equations using the Taylor series method:
x ² + y ¢ - 4 x + 12 = 0
y ² - 10 x ¢ - y + 7 = 0
and
x ² + 2 x ¢ + 2 y ¢ + 3z ¢ + x = 1
y¢ + z¢ - x = 0
x ¢ + z = 0.
205
www.pdfgrip.com
Chapter 8 ■ Differential Equations
>> pretty(simple(sym(maple('dsolve({diff(x(t),t$2)+diff(y(t),t)-4*x+12=0,
diff(y(t),t$2)-10*diff(x(t),t)-y(t)+7=0,x(0)=1,y(0)=1,D(x)(0)=1,D(y)(0)=1},
{x(t),y(t)},series)'))))
2 3 4 147 5 6
{y(t) = 1 + t + 2 t - 89/6 t + 1/6 t + --- t + O(t ),
40
2 53 4 5 6
x(t) = 1 + t - 9/2 t + -- t - 1/30 t + O(t )}
24
>> pretty(simple(sym(maple('dsolve({diff(x(t),t$2)+2*diff(x(t),t)+2*diff(y(t),t)+ +3*diff
(z(t),t)+x(t)=1,diff(y(t),t)+diff(z(t),t)-x(t)=0,diff(x(t),t)+z(t)=0}, {x(t),y(t),z(t)},
series)'))))
2
{x(t) = x(0) + D(x)(0) t + (- D(x)(0) - 1/2 x(0) + 1/2) t
3 4
+ (1/2 D(x)(0) + 1/3 x(0) - 1/3) t + (- 1/6 D(x)(0) - 1/8 x(0) + 1/8) t
5 6
+ (1/24 D(x)(0) + 1/30 x(0) - 1/30) t + O(t ),
2 3
z(t) = z(0) + x(0) t + 1/2 D(x)(0) t + (- 1/3 D(x)(0) - 1/6 x(0) + 1/6) t
4
+ (1/8 D(x)(0) + 1/12 x(0) - 1/12) t
5 6
+ (- 1/30 D(x)(0) - 1/40 x(0) + 1/40) t + O(t ),
2 3
y(t) = y(0) + x(0) t + 1/2 D(x)(0) t + (- 1/3 D(x)(0) - 1/6 x(0) + 1/6) t
4
+ (1/8 D(x)(0) + 1/12 x(0) - 1/12) t
5 6
+ (- 1/30 D(x)(0) - 1/40 x(0) + 1/40) t + O(t ) }
206
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-15
Figure 8-4.
207
www.pdfgrip.com
Chapter 8 ■ Differential Equations
We find the degree 2 polynomial which is the best fit to the set of solution points. The equation of this parabola
will be an approximate solution of the differential equation.
>> pretty(vpa(poly2sym(polyfit((-0.3:.1:0.3),y,2))))
2
.5747427827483573 x + 1.041293962469090 x +.4991457846921903
This yields a degree 2 polynomial approximation to the solution y(x) of the equation.
208
www.pdfgrip.com
Chapter 8 ■ Differential Equations
209
www.pdfgrip.com
Chapter 8 ■ Differential Equations
constcoeffsol(reqn,fnc(var),{cond1,...,condm},output=basis)
Returns a basis for the solutions of reqn
hypergeomsols(reqn,fnc(var),{cond1,...,condm})
Returns hypergeometric solutions of reqn. Also supports the option output =
basis. With the option output = onesol gives a simple solution.
polysols (reqn, fnc (var), {cond1,..., condm})
Returns polynomial solutions of reqn. Supports the additional options
output = basis and output = onesol.
ratpolysols (reqn, fnc (var), {cond1,..., condm})
Returns rational polynomial solutions of reqn with the given initial conditions
REcreate(reqn,fnc(var),{cond1,...,condn})
Creates a recurrence equation for reqn
REcreate ({reqn1,..., reqnn}, {(var) fnc1,..., fncn (var)}, {cond1,..., condm})
Creates n recurrence equations
REcontent (reqn, fnc (var), {cond1,..., condm})
Returns the content of the recurrence equation reqn
REprimpart(reqn,fnc(var),{cond1,...,condm})
Returns the primitive part of the recurrence equation
REreduceorder (reqn, fnc (var), {cond1,..., condm}, expr)
Reduction of order for the recurrent equation reqn where expr is a partial solution
REreduceorder (reqn, fnc (var), {cond1,..., condm}, [expr1,..., exprn])
Specifies n partial solutions
REtoDE (reqn, fnc (var), {cond1,..., condm})
Converts a recurrence equation reqn to a differential equation
REtodelta (reqn, fnc (var), {cond1,..., condm})
Converts the recurrence equation reqn to a difference equation
REtoproc(reqn,fnc(var),{cond1,...,condm})
Converts the recurrence equation reqn into a procedure
shift(expr, var, int)
Coincides with subs(var=var+int,expr)
shift(expr, var)
Takes int = 1 in the above by default
210
www.pdfgrip.com
Chapter 8 ■ Differential Equations
211
www.pdfgrip.com
Chapter 8 ■ Differential Equations
Exercise 8-16
Exercise 8-17
Find the general term of the sequences of real numbers defined by the following recurrence laws:
xn −nxm* xn *xn + 1 = xn + 1 x0 = 1,
yn + 2 −2yn + 1 + 5yn = cos(3n) y0 = y1 = 1
>> pretty(sym(maple('rsolve({x(n)-n*x(n)*x(n+1)=x(n+1),x(0)=1},x) ')))
2
----------
2
n - n + 2
>> pretty(sym(maple('simplify(evalf(rsolve({x(n+2)-2*x(n+1)+5*x(n)=
=cos(3*n), x(0)=1,x(1)=1},x))) ')))
.4373424525 exp((.8047189562 - 1.107148718 I) n)
+ .4373424525 exp((.8047189562 + 1.107148718 I) n)
212
www.pdfgrip.com
Chapter 8 ■ Differential Equations
3 3
- .9160727930 cos(n - 3.) + .4362567143 cos(n - 5.)
3 3
+ .1725076052 cos(n - 2.) - .1986350636 cos(n - 4.)
- .06265675740 I exp((.8047189562 + 1.107148718 I) n)
+ .06265675750 I exp((.8047189562 - 1.107148718 I) n)
+ .7931668115 cos(n - 3.) + .07520876265 cos(n - 2.)
213
www.pdfgrip.com
Contents
v
www.pdfgrip.com
■ Contents
vii
www.pdfgrip.com
■ Contents
viii
www.pdfgrip.com
César Pérez López is a Professor at the Department of Statistics and Operations Research at the University of
Madrid. César is also a Mathematician and Economist at the National Statistics Institute (INE) in Madrid, a body
which belongs to the Superior Systems and Information Technology Department of the Spanish Government.
César also currently works at the Institute for Fiscal Studies in Madrid.
ix
www.pdfgrip.com
Coming Soon
xi