Calculus

Differentiation

To illustrate how to take derivatives using Symbolic Math Toolbox software, first create a symbolic expression:

syms x
f = sin(5*x)

The command

diff(f)

differentiates f with respect to x:

ans =
5*cos(5*x)

As another example, let

g = exp(x)*cos(x)

where exp(x) denotes e^x, and differentiate g:

diff(g)
ans =
exp(x)*cos(x) - exp(x)*sin(x)

To take the second derivative of g, enter

diff(g,2)
ans =
(-2)*exp(x)*sin(x)

You can get the same result by taking the derivative twice:

diff(diff(g))
ans =
(-2)*exp(x)*sin(x)

In this example, MATLAB software automatically simplifies the answer. However, in some cases, MATLAB might not simply an answer, in which case you can use the simplify command. For an example of this, see More Examples.

Note that to take the derivative of a constant, you must first define the constant as a symbolic expression. For example, entering

c = sym('5');
diff(c)

returns

ans =
0

If you just enter

diff(5)

MATLAB returns

ans =
     []

because 5 is not a symbolic expression.

Derivatives of Expressions with Several Variables

To differentiate an expression that contains more than one symbolic variable, you must specify the variable that you want to differentiate with respect to. The diff command then calculates the partial derivative of the expression with respect to that variable. For example, given the symbolic expression

syms s t
f = sin(s*t)

the command

diff(f,t)

calculates the partial derivative . The result is

ans = 
s*cos(s*t)

To differentiate f with respect to the variable s, enter

diff(f,s)

which returns:

ans =
t*cos(s*t)

If you do not specify a variable to differentiate with respect to, MATLAB chooses a default variable. Basically, the default variable is the letter closest to x in the alphabet. See the complete set of rules in Finding a Default Symbolic Variable. In the preceding example, diff(f) takes the derivative of f with respect to t because the letter t is closer to x in the alphabet than the letter s is. To determine the default variable that MATLAB differentiates with respect to, use the symvar command:

symvar(f, 1)

ans = 
t

To calculate the second derivative of f with respect to t, enter

diff(f, t, 2)

which returns

ans =
-s^2*sin(s*t)

Note that diff(f, 2) returns the same answer because t is the default variable.

More Examples

To further illustrate the diff command, define a, b, x, n, t, and theta in the MATLAB workspace by entering

syms a b x n t theta

The table below illustrates the results of entering diff(f).

syms x n;
f = x^n;

diff(f)

ans =
n*x^(n - 1)

syms a b t;
f = sin(a*t + b);

diff(f)

ans =
a*cos(b + a*t)

syms theta;
f = exp(i*theta);

diff(f)

ans =
i*exp(i*theta)

To differentiate the Bessel function of the first kind,besselj(nu,z), with respect to z, type

syms nu z
b = besselj(nu,z);
db = diff(b)

which returns

db =
(nu*besselj(nu, z))/z - besselj(nu + 1, z)

The diff function can also take a symbolic matrix as its input. In this case, the differentiation is done element-by-element. Consider the example

syms a x
A = [cos(a*x),sin(a*x);-sin(a*x),cos(a*x)]

which returns

A =
[  cos(a*x), sin(a*x)]
[ -sin(a*x), cos(a*x)]

The command

diff(A)

returns

ans =
[ -a*sin(a*x),  a*cos(a*x)]
[ -a*cos(a*x), -a*sin(a*x)]

You can also perform differentiation of a vector function with respect to a vector argument. Consider the transformation from Euclidean (x, y, z) to spherical coordinates as given by , , and . Note that corresponds to elevation or latitude while denotes azimuth or longitude.

To calculate the Jacobian matrix, J, of this transformation, use the jacobian function. The mathematical notation for J is

For the purposes of toolbox syntax, use l for and f for . The commands

syms r l f
x = r*cos(l)*cos(f); y = r*cos(l)*sin(f); z = r*sin(l);
J = jacobian([x; y; z], [r l f])

return the Jacobian

J =
[ cos(f)*cos(l), -r*cos(f)*sin(l), -r*cos(l)*sin(f)]
[ cos(l)*sin(f), -r*sin(f)*sin(l),  r*cos(f)*cos(l)]
[        sin(l),         r*cos(l),                0]

and the command

detJ = simple(det(J))

returns

detJ = 
-r^2*cos(l)

The arguments of the jacobian function can be column or row vectors. Moreover, since the determinant of the Jacobian is a rather complicated trigonometric expression, you can use the simple command to make trigonometric substitutions and reductions (simplifications). The section Simplifications and Substitutions discusses simplification in more detail.

A table summarizing diff and jacobian follows.

Back to Top

Limits

The fundamental idea in calculus is to make calculations on functions as a variable "gets close to" or approaches a certain value. Recall that the definition of the derivative is given by a limit

provided this limit exists. Symbolic Math Toolbox software enables you to calculate the limits of functions directly. The commands

syms h n x
limit((cos(x+h) - cos(x))/h, h, 0)

which return

ans =
-sin(x)

and

limit((1 + x/n)^n, n, inf)

which returns

ans =
exp(x)

illustrate two of the most important limits in mathematics: the derivative (in this case of cos(x)) and the exponential function.

One-Sided Limits

You can also calculate one-sided limits with Symbolic Math Toolbox software. For example, you can calculate the limit of x/|x|, whose graph is shown in the following figure, as x approaches 0 from the left or from the right.

To calculate the limit as x approaches 0 from the left,

enter

syms x;
limit(x/abs(x), x, 0, 'left')

This returns

ans =
 -1

To calculate the limit as x approaches 0 from the right,

enter

syms x;
limit(x/abs(x), x, 0, 'right')

This returns

ans =
1

Since the limit from the left does not equal the limit from the right, the two- sided limit does not exist. In the case of undefined limits, MATLAB returns NaN (not a number). For example,

syms x;
limit(x/abs(x), x, 0)

returns

ans =
NaN

Observe that the default case, limit(f) is the same as limit(f,x,0). Explore the options for the limit command in this table, where f is a function of the symbolic object x.

Back to Top

Integration

If f is a symbolic expression, then

int(f)

attempts to find another symbolic expression, F, so that diff(F) = f. That is, int(f) returns the indefinite integral or antiderivative of f (provided one exists in closed form). Similar to differentiation,

int(f,v)

uses the symbolic object v as the variable of integration, rather than the variable determined by symvar. See how int works by looking at this table.

In contrast to differentiation, symbolic integration is a more complicated task. A number of difficulties can arise in computing the integral:

• The antiderivative, F, may not exist in closed form.

• The antiderivative may define an unfamiliar function.

• The antiderivative may exist, but the software can't find it.

• The software could find the antiderivative on a larger computer, but runs out of time or memory on the available machine.

Nevertheless, in many cases, MATLAB can perform symbolic integration successfully. For example, create the symbolic variables

syms a b theta x y n u z

The following table illustrates integration of expressions containing those variables.

syms x n;
f = x^n;

int(f)

ans =
piecewise([n = -1, log(x)], [n <> -1, x^(n + 1)/(n + 1)])

syms y;
f = y^(-1);

int(f)

ans =
log(y)

syms x n;
f = n^x;

int(f)

ans =
n^x/log(n)

syms a b theta;
f = sin(a*theta+b);

int(f)

ans =
-cos(b + a*theta)/a

syms u;
f = 1/(1+u^2);

int(f)

ans =
atan(u)

syms x;
f = exp(-x^2);

int(f)

ans =
(pi^(1/2)*erf(x))/2

In the last example, exp(-x^2), there is no formula for the integral involving standard calculus expressions, such as trigonometric and exponential functions. In this case, MATLAB returns an answer in terms of the error function erf.

If MATLAB is unable to find an answer to the integral of a function f, it just returns int(f).

Definite integration is also possible.

Here are some additional examples.

syms x;
f = x^7;

a = 0;
b = 1;

int(f, a, b)

ans =
1/8

syms x;
f = 1/x;

a = 1;
b = 2;

int(f, a, b)

ans =
log(2)

syms x;
f = log(x)*sqrt(x);

a = 0;
b = 1;

int(f, a, b)

ans =
-4/9

syms x;
f = exp(-x^2);

a = 0;
b = inf;

int(f, a, b)

ans =
pi^(1/2)/2

syms z;
f = besselj(1,z)^2;

a = 0;
b = 1;

int(f, a, b)

ans =
hypergeom([3/2, 3/2], [2, 5/2, 3], -1)/12

For the Bessel function (besselj) example, it is possible to compute a numerical approximation to the value of the integral, using the double function. The commands

syms z
a = int(besselj(1,z)^2,0,1)

return

a =
hypergeom([3/2, 3/2], [2, 5/2, 3], -1)/12

and the command

a = double(a)

returns

a =
    0.0717

Integration with Real Parameters

One of the subtleties involved in symbolic integration is the "value" of various parameters. For example, if a is any positive real number, the expression

is the positive, bell shaped curve that tends to 0 as x tends to ±∞. You can create an example of this curve, for a = 1/2, using the following commands:

syms x
a = sym(1/2);
f = exp(-a*x^2);
ezplot(f)

However, if you try to calculate the integral

without assigning a value to a, MATLAB assumes that a represents a complex number, and therefore returns a piecewise answer that depends on the argument of a. If you are only interested in the case when a is a positive real number, you can calculate the integral as follows:

syms a positive;

The argument positive in the syms command restricts a to have positive values. Now you can calculate the preceding integral using the commands

syms x;
f = exp(-a*x^2);
int(f, x, -inf, inf)

This returns

ans =
pi^(1/2)/a^(1/2)

Integration with Complex Parameters

To calculate the integral

for complex values of a, enter

syms a x clear 
f = 1/(a^2 + x^2);
F = int(f, x, -inf, inf)

syms is used with the clear option to clear the real property that was assigned to a in the preceding example — see Deleting Symbolic Objects and Their Assumptions.

The preceding commands produce the complex output

F = 
(pi*signIm(i/a))/a

The function signIm is defined as:

To evaluate F at a = 1 + i, enter

g = subs(F, 1 + i)

g = 
pi/(2*i)^(1/2)

double(g)

ans =
   1.5708 - 1.5708i

Back to Top

Symbolic Summation

You can compute symbolic summations, when they exist, by using the symsum command. For example, the p-series

sums to , while the geometric series

1 + x + x² + ...

sums to 1/(1 – x), provided . These summations are demonstrated below:

syms x k
s1 = symsum(1/k^2, 1, inf)
s2 = symsum(x^k, k, 0, inf)

s1 =
pi^2/6

s2 =
piecewise([1 <= x, Inf], [abs(x) < 1, -1/(x - 1)])

Back to Top

Taylor Series

The statements

syms x
f = 1/(5 + 4*cos(x));
T = taylor(f, 8)

return

T =
(49*x^6)/131220 + (5*x^4)/1458 + (2*x^2)/81 + 1/9

which is all the terms up to, but not including, order eight in the Taylor series for f(x):

Technically, T is a Maclaurin series, since its base point is a = 0.

The command

pretty(T)

prints T in a format resembling typeset mathematics:

 
      6       4      2 
  49 x     5 x    2 x    1 
  ------ + ---- + ---- + - 
  131220   1458    81    9

These commands

syms x
g = exp(x*sin(x))
t = taylor(g, 12, 2);

generate the first 12 nonzero terms of the Taylor series for g about x = 2.

t is a large expression; enter

size(char(t))

ans =
           1       109959

to find that t has more than 100,000 characters in its printed form. In order to proceed with using t, first simplify its presentation:

t = simplify(t);
size(char(t))

ans =
           1       11585

To simplify t even further, use the simple function:

t = simple(t);
size(char(t))

ans =
           1        6988

Next, plot these functions together to see how well this Taylor approximation compares to the actual function g:

xd = 1:0.05:3; yd = subs(g,x,xd);
ezplot(t, [1, 3]); hold on;
plot(xd, yd, 'r-.')
title('Taylor approximation vs. actual function');
legend('Taylor','Function')

Special thanks to Professor Gunnar Bäckstrøm of UMEA in Sweden for this example.

Back to Top

Calculus Example

This section describes how to analyze a simple function to find its asymptotes, maximum, minimum, and inflection point. The section covers the following topics:

• Defining the Function

• Finding the Asymptotes

• Finding the Maximum and Minimum

• Finding the Inflection Point

Defining the Function

The function in this example is

To create the function, enter the following commands:

syms x
num = 3*x^2 + 6*x -1;
denom = x^2 + x - 3;
f = num/denom

This returns

f = 
(3*x^2 + 6*x - 1)/(x^2 + x - 3)

You can plot the graph of f by entering

ezplot(f)

This displays the following plot.

Finding the Asymptotes

To find the horizontal asymptote of the graph of f, take the limit of f as x approaches positive infinity:

limit(f, inf)
ans = 
3

The limit as x approaches negative infinity is also 3. This tells you that the line y = 3 is a horizontal asymptote to the graph.

To find the vertical asymptotes of f, set the denominator equal to 0 and solve by entering the following command:

roots = solve(denom)

This returns to solutions to :

roots =
 
   13^(1/2)/2 - 1/2
 - 13^(1/2)/2 - 1/2

This tells you that vertical asymptotes are the lines

and

You can plot the horizontal and vertical asymptotes with the following commands:

ezplot(f)
hold on % Keep the graph of f in the figure
% Plot horizontal asymptote
plot([-2*pi 2*pi], [3 3],'g')
% Plot vertical asymptotes
plot(double(roots(1))*[1 1], [-5 10],'r')
plot(double(roots(2))*[1 1], [-5 10],'r')
title('Horizontal and Vertical Asymptotes')
hold off

Note that roots must be converted to double to use the plot command.

The preceding commands display the following figure.

To recover the graph of f without the asymptotes, enter

ezplot(f)

Finding the Maximum and Minimum

You can see from the graph that f has a local maximum somewhere between the points x = –2 and x = 0, and might have a local minimum between x = –6 and x = –2. To find the x-coordinates of the maximum and minimum, first take the derivative of f:

f1 = diff(f)

This returns

f1 = 
(6*x + 6)/(x^2 + x - 3) - ((2*x + 1)*(3*x^2 + 6*x - 1))/(x^2 + x - 3)^2

To simplify this expression, enter

f1 = simplify(f1)

which returns

f1 =
 -(3*x^2 + 16*x + 17)/(x^2 + x - 3)^2

You can display f1 in a more readable form by entering

pretty(f1)

which returns

     
       2 
    3 x  + 16 x + 17 
  - ---------------- 
       2         2 
     (x  + x - 3)

Next, set the derivative equal to 0 and solve for the critical points:

crit_pts = solve(f1)

This returns

crit_pts =
 
   13^(1/2)/3 - 8/3
 - 13^(1/2)/3 - 8/3

It is clear from the graph of f that it has a local minimum at

and a local maximum at

Note   MATLAB does not always return the roots to an equation in the same order.

You can plot the maximum and minimum of f with the following commands:

ezplot(f)
hold on
plot(double(crit_pts), double(subs(f,crit_pts)),'ro')
title('Maximum and Minimum of f')
text(-5.5,3.2,'Local minimum')
text(-2.5,2,'Local maximum')
hold off

This displays the following figure.

Finding the Inflection Point

To find the inflection point of f, set the second derivative equal to 0 and solve.

f2 = diff(f1);
inflec_pt = solve(f2);
double(inflec_pt)

This returns

ans =
  -5.2635          
  -1.3682 - 0.8511i
  -1.3682 + 0.8511i

In this example, only the first entry is a real number, so this is the only inflection point. (Note that in other examples, the real solutions might not be the first entries of the answer.) Since you are only interested in the real solutions, you can discard the last two entries, which are complex numbers.

inflec_pt = inflec_pt(1)

To see the symbolic expression for the inflection point, enter

pretty(simplify(inflec_pt))

This returns

 
                           /              1/2 \1/3 
             13            |          2197    |      8 
  - -------------------- - | 169/54 - ------- |    - - 
                       1   \             18   /      3 
                       - 
      /           1/2 \3 
      | 169   2197    | 
    9 | --- - ------- | 
      \ 54      18    /

To plot the inflection point, enter

ezplot(f, [-9 6])
hold on
plot(double(inflec_pt), double(subs(f,inflec_pt)),'ro')
title('Inflection Point of f')
text(-7,2,'Inflection point')
hold off

The extra argument, [-9 6], in ezplot extends the range of x values in the plot so that you see the inflection point more clearly, as shown in the following figure.

Back to Top

Extended Calculus Example

This section presents an extended example that illustrates how to find the maxima and minima of a function. The section covers the following topics:

• Defining the Function

• Finding the Zeros of f3

• Finding the Maxima and Minima of f2

• Integrating

Defining the Function

The starting point for the example is the function

You can create the function with the commands

syms x 
f = 1/(5+4*cos(x))

which return

f =
1/(4*cos(x) + 5)

The example shows how to find the maximum and minimum of the second derivative of f(x). To compute the second derivative, enter

f2 = diff(f, 2)

which returns

f2 =
(4*cos(x))/(4*cos(x) + 5)^2 + (32*sin(x)^2)/(4*cos(x) + 5)^3

Equivalently, you can type f2 = diff(f, x, 2). The default scaling in ezplot cuts off part of the graph of f2. You can set the axes limits manually to see the entire function:

ezplot(f2)   
axis([-2*pi 2*pi -5 2])
title('Graph of f2')

From the graph, it appears that the maximum value of is 1 and the minimum value is -4. As you will see, this is not quite true. To find the exact values of the maximum and minimum, you only need to find the maximum and minimum on the interval (–π, π]. This is true because is periodic with period 2π, so that the maxima and minima are simply repeated in each translation of this interval by an integer multiple of 2π. The next two sections explain how to do find the maxima and minima.

Finding the Zeros of f3

The maxima and minima of occur at the zeros of . The statements

f3 = diff(f2);
pretty(f3)

compute and display it in a more readable form:

 
              3 
    384 sin(x)         4 sin(x)       96 cos(x) sin(x) 
  --------------- - --------------- + ---------------- 
                4                 2                 3 
  (4 cos(x) + 5)    (4 cos(x) + 5)    (4 cos(x) + 5)

You can simplify this expression using the statements

f3 = simple(f3);
pretty(f3)

 
                       2 
  4 sin(x) (- 16 cos(x)  + 80 cos(x) + 71) 
  ---------------------------------------- 
                            4 
              (4 cos(x) + 5)

Now, to find the zeros of , enter

zeros = solve(f3)

This returns a 5-by-1 symbolic matrix

zeros = 
 acos(5/2 - (3*19^(1/2))/4)
 acos((3*19^(1/2))/4 + 5/2)
                          0
  -acos(5/2 - 3/4*19^(1/2))
  -acos(3/4*19^(1/2) + 5/2)

each of whose entries is a zero of . The commands

format;
% Default format of 5 digits
zerosd = double(zeros)

convert the zeros to double form:

zerosd =
   2.4483          
        0 + 2.4381i
        0          
  -2.4483          
        0 - 2.4381i        

So far, you have found three real zeros and two complex zeros. However, as the following graph of f3 shows, these are not all its zeros:

ezplot(f3)
hold on;
plot(zerosd,0*zerosd,'ro') % Plot zeros
plot([-2*pi,2*pi], [0,0],'g-.'); % Plot x-axis
title('Graph of f3')

The red circles in the graph correspond to zerosd(1), zerosd(3), and zerosd(4). As you can see in the graph, there are also zeros at ±π. The additional zeros occur because contains a factor of sin(x), which is zero at integer multiples of π. The function, solve(sin(x)), however, only finds the zero at x = 0.

A complete list of the zeros of in the interval (–π, π] is

zerosd = [zerosd(1) zerosd(3) zerosd(4) pi];

You can display these zeros on the graph of with the following commands:

ezplot(f3)
hold on;
plot(zerosd,0*zerosd,'ro')
plot([-2*pi,2*pi], [0,0],'g-.'); % Plot x-axis
title('Zeros of f3')
hold off;

Finding the Maxima and Minima of f2

To find the maxima and minima of , calculate the value of at each of the zeros of . To do so, substitute zeros into f2 and display the result below zeros:

[zerosd; subs(f2,zerosd)]

ans =
    2.4483         0   -2.4483    3.1416
    1.0051    0.0494    1.0051   -4.0000

This shows the following:

• has an absolute maximum at x = ±2.4483, whose value is 1.0051.

• has an absolute minimum at x = π, whose value is -4.

• has a local minimum at x = 0, whose value is 0.0494.

You can display the maxima and minima with the following commands:

clf
ezplot(f2)
axis([-2*pi 2*pi -4.5 1.5])
ylabel('f2');
title('Maxima and Minima of f2')
hold on
plot(zerosd, subs(f2,zerosd), 'ro')
text(-4, 1.25, 'Absolute maximum')
text(-1,-0.25,'Local minimum')
text(.9, 1.25, 'Absolute maximum')
text(1.6, -4.25, 'Absolute minimum')
hold off;

This displays the following figure.

The preceding analysis shows that the actual range of is [–4, 1.0051].

Integrating

Integrate f(x):

F = int(f)

The result

F =
(2*atan(tan(x/2)/3))/3

involves the arctangent function.

Note that F(x) is not an antiderivative of f(x) for all real numbers, since it is discontinuous at odd multiples of π, where tan (x) is singular. You can see the gaps in F(x) in the following figure.

ezplot(F)

To change F(x) into a true antiderivative of f(x) that is differentiable everywhere, you can add a step function to F(x). The height of the steps is the height of the gaps in the graph of F(x). You can determine the height of the gaps by taking the limits of F(x) as x approaches π from the left and from the right. The limit from the left is

limit(F, x, pi, 'left')

ans =
pi/3

On the other hand, the limit from the right is

limit(F, x, pi, 'right')

ans =
-pi/3

The height of the gap is the distance between the left-and right-hand limits, which is 2π/3 as shown in the following figure.

You can create the step function using the round function, which rounds numbers to the nearest integer, as follows:

J = sym(2*pi/3)*sym('round(x/(2*pi))');

Each step has width 2π, and the jump from one step to the next is 2π/3, as shown in the following figure, generated with

ezplot(J, [-2*pi 2*pi])

Next, add the step function J(x) to F(x) with the following code:

F1 = F + J

F1 =
(2*atan(tan(x/2)/3))/3 + (2*pi*round(x/(2*pi)))/3

Adding the step function raises the section of the graph of F(x) on the interval [π, 3π) up by 2π/3, lowers the section on the interval (–3π, –π] down by 2π/3, and so on, as shown in the following figure.

When you plot the result by entering

ezplot(F1)

you see that this representation does have a continuous graph.

Back to Top

Using Symbolic Math Toolbox Software

Simplifications and Substitutions

Recommended Products

• MATLAB
• Signal Processing Toolbox
• Optimization Toolbox
• Statistics Toolbox
• Control System Toolbox

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit