5. Using the Symbolic Math Toolbox¶

The majority of our consideration relates to numerical algorithms. But here, we focus on analytic solutions to mathematical problems. We seek analytical solutions in mathematics courses, and we typically solve them with pencil and paper. Computers are, of course, more suited to numerical algorithms. But MATLAB does have the ability to find analytic solutions when they are available. When a computer solves a mathematics problem analytically, we say that the computer is operating on symbols in addition to numbers. With symbolic operations, the computer would simplify expressions such as $6/3$ to $2$ , but would leave numbers like $\pi$ , $e$ , and $17/7$ as irreducible numbers.

The Symbolic MATH Toolbox (SMT) is required to access the features of MATLAB discussed here. The SMT is part of the student MATLAB bundle, but is an additional purchase with other MATLAB licenses.

5.1. Throwing a Ball Up¶

As a quick tour of some of the features of the SMT, consider throwing a ball straight upward. The scene is depicted in Fig. 5.1. The algebraic equations describing the force of gravity on a moving ball are readily available on the web or in physics and dynamics textbooks. But let us use the SMT to derive the equations from differential equations and apply them. In this example, we will use functions from the SMT with minimal comments about them. Then we will give more details about the SMT functions beginning in Symbolic Algebra.

Fig. 5.1 A ball is thrown upward from height $y_0$ with an initial velocity of $V_0$ . When the ball reaches its maximum height, $y_{max}$ at time $t_{max}$ , its velocity is zero, $v(t_{max}) = 0$ . The ball hits the ground at time $t_{final}$ , $y(t_{final}) = 0$ .

The vertical position of the ball with respect to time is $y(t)$ . The velocity of the ball is the derivative of the position,

(5.1)¶ $v(t) = \frac{dy(t)}{dt}.$

The derivative of the velocity is the acceleration, which in this case is a constant from gravity in the negative direction. The acceleration is the second derivative of the position.

(5.2)¶ $a(t) = \frac{dv(t)}{dt} = \frac{dy^2(t)}{dt^2} = g$

$g = -9.8 \; m/sec^2 \mbox{ or } g = -32 \; feet/sec^2$

Our starting point, equation (5.2), is a second order differential equation. In addition to describing the differential equation, we need to describe the initial conditions of the problem. The initial height of the ball is given by $y(0) = y_0$ ; and the initial velocity from the throw is $v(0) = V_0$ .

Symbolic variables and functions are first created in MATLAB with the syms statement. Another way to create a symbolic variable is to pass a value or expression to the sym function, such as z = sym(5). The sym function returns a symbolic variable.

The differential equations in equation (5.1) and equation (5.2) are solved in the following MATLAB code. We will declare y(t) (position) as a symbolic function, and y0 (initial position), V0 (initial velocity), and g (gravity) as symbolic variables. Time, t, is implicitly declared from y(t). MATLAB solves symbolic differential equations with the dsolve function.

Note in the code below that the equal sign (=) is, as always, used for assignment statements, but the double equal sign (==) is used to declare equality relationships in equations.

>> syms y(t) y0 V0 g
>> diffeq = diff(y(t),t,2) == g
diffeq =
    diff(y(t), t, 2) == g
>> v(t) = diff(y(t),t)
v(t) =
    diff(y(t), t)

>> cond = [y(0) == y0, v(0) == V0]
cond =
    [ y(0) == y0, subs(diff(y(t), t), t, 0) == V0]
>> y(t) = dsolve(diffeq, cond)
y(t) =
    (g*t^2)/2 + V0*t + y0

The diff function takes the derivative of the position to find the velocity.

>> v(t) = diff(y(t),t)
v(t) =
    V0 + g*t

The maximum height of the ball occurs when the velocity is zero. The solve function finds solutions to algebraic equations.

>> t_max = solve(v(t) == 0)
t_max =
    -V0/g
>> y_max = y(t_max)
y_max =
    -V0^2/(2*g) + y0

The ball lands on the ground when y(t) is zero. Two answers are returned because y(t) is a quadratic equation.

>> t_final = solve(y(t) == 0, t)
t_final =
    -(V0 + (V0^2 - 2*g*y0)^(1/2))/g
    -(V0 - (V0^2 - 2*g*y0)^(1/2))/g

We can give the variables numerical values and substitute (subs) the values into the equations, and see a decimal version of the number with vpa.

% From 1 meter off groud, throw up at 10 m/s
>> y0 = 1; V0 = 10; g = -9.8;
>> subs(t_max)             % time to max height
ans =
    50/49
>> vpa(subs(t_max))
ans =
    1.0204081632653061224489795918367 % seconds

>> subs(y_max)
ans =
    299/49
>> vpa(subs(y_max))
ans =
    6.1020408163265306122448979591837 % meters

>> subs(t_final(1))           % time until the ball hits the ground
ans =
    (5^(1/2)*598^(1/2))/49 + 50/49
>> vpa(subs(t_final(1)))
ans =
    2.1363447440401890728050707611418 % seconds

5.2. Symbolic Algebra¶

5.2.1. Collect¶

The collect function combines variables of the same power of the specified variable in an expression.

>> syms x y
>> f = (2*x^2 + 3*y) * (x^3 + x^2 + x*y^2 + 2*y);
>> collect(f,x)
ans =
    2*x^5 + 2*x^4 + (2*y^2 + 3*y)*x^3 + 7*y*x^2 + 3*y^3*x + 6*y^2
>> collect(f,y)
ans =
    3*x*y^3 + (2*x^3 + 6)*y^2 + (3*x^3 + 7*x^2)*y + 2*x^2*(x^3 + x^2)

5.2.2. Factor¶

The factor function finds the set of irreducible polynomials whose product form the higher order polynomial argument to factor. The factoring operation is useful for finding the roots of the polynomial.

>> g = x^4 + 3*x^3 - 19*x^2 - 27*x + 90;
>> factor(g)
ans =
    [ x - 2, x - 3, x + 5, x + 3]

5.2.3. Expand¶

The expand function multiplies polynomial products to get a single polynomial.

>> f = (3*x^2 + 2*x)*(2*x - 5);
>> expand(f)
ans =
    6*x^3 - 11*x^2 - 10*x
>> g
g =
    (x - 2)*(x - 3)*(x + 3)*(x + 5)
>> expand(g)
ans =
    x^4 + 3*x^3 - 19*x^2 - 27*x + 90

5.2.4. Simplify¶

The simplify function generates a simpler form of an expression.

>> g
g =
    x^4 + 3*x^3 - 19*x^2 - 27*x + 90
>> h = expand((x-2)*(x+5))
h =
    x^2 + 3*x - 10
>> k = g/h
k =
    (x^4 + 3*x^3 - 19*x^2 - 27*x + 90)/(x^2 + 3*x - 10)
>> simplify(k)
ans =
    x^2 - 9

>> f = cos(x)*cos(y) + sin(x)*sin(y);
>> simplify(f)
ans =
    cos(x - y)

>> f = cos(x)^2 + sin(x)^2;
>> simplify(f)
ans =
    1

5.2.5. Solve¶

Finding the roots of equations, i.e., $x$ , where $f(x) = 0$ , is addressed here and also in the context of two other topics.

To find the numerical roots of a polynomial using a linear algebra algorithm, see the roots function described in Roots of a Polynomial by Eigenvalues.
To use numerical methods to find the roots of any real equation, including non-polynomials, see the fzero function described in Numerical Roots of a Function.

The SMT function finds the roots of any function. The solve(s,x) function returns the set of symbolic values of x that satisfy the expression s. The default is to solve for s == 0, but an equality statement may be used as part of the expression. Be certain to use an equality statement (==), not an assignment statement (=).

>> g(x) = x^2 - 1;
>> solve(g(x),x)
ans =
    -1
     1

>> solve(g(x) == 4,x)
ans =
     5^(1/2)
    -5^(1/2)

>> h(x) = x + 2;
>> solve(g(x) == h(x), x)
ans =
    1/2 - 13^(1/2)/2
    13^(1/2)/2 + 1/2

>> syms t y(t)
>> y(t) = cos(pi*t)  % non-polynomial
y(t) =
    cos(pi*t)
>> solve(y(t), t)
ans =
    1/2

>> y(t) = 1 + t^2;
>> solve(y(t), t)
ans =              % even complex roots
   -1i
    1i

5.2.6. Subs¶

The subs(s) function returns a copy of s after replacing the symbolic variables in s with their values from the MATLAB workspace. The subs function was demonstrated in Throwing a Ball Up.

>> syms g(x) a b x
>> g(x) = a*x^2 + b*x;
>> a = 2; b = 6;
>> subs(g)
ans(x) =
    2*x^2 + 6*x
>> subs(g(3))
ans =
    36
>> syms c
>> h = g/c
h(x) =
    (a*x^2 + b*x)/c
>> c = 3;
>> subs(h(3))
ans =
    12

5.2.7. Vpa¶

The Symbolic Math Toolbox tries to keep variables as exact symbolic values rather than as decimal approximations. The vpa function converts symbolic expressions to decimal values. The vpa function gives decimal answers. The name vpa is an acronym for variable–precision floating-point arithmetic. As the name implies, the vpa(x, d) function evaluates x to at least d digits. The default number of digits is 32. The maximum number of digits is $2^{29}$ .

>> g(x) = x^2 - 1;
>> h = solve(g(x) == 4, x)
h =
    5^(1/2)
   -5^(1/2)
>> vpa(h)
ans =
    2.2360679774997896964091736687313
   -2.2360679774997896964091736687313

5.3. Symbolic Calculus¶

5.3.1. Symbolic Derivatives¶

The diff(s, var, n) function takes the symbolic derivative of an expression with respect to a specified variable. An optional third argument specifies the order of the derivative. As with many SMT functions, the default variable is x. The variable is needed when the equation has more than one variable, but it never hurts to specify it.

>> syms x
>> diff(x^2 + 3*x + 2,x)
ans =
    2*x + 3
>> diff(x^2 + 3*x + 2,x,2)  % 2nd derivative
ans =
    2
>> diff(x*sin(3*x))
ans =
    sin(3*x) + 3*x*cos(3*x)
>> diff(exp(-x/2))
ans =
   -exp(-x/2)/2

5.3.2. Symbolic Integration¶

The int function performs integration of an expression. The variable specification is optional if the expression has only one symbolic variable. Thus valid forms for indefinite integrals are int(s) and int(s,var).

$\int x^2 + 3\,x + 2 \, dx = \frac{1}{3} x^3 + \frac{3}{2} x^2 + 2\,x$

>> syms x t
>> int(x^2 + 3*x + 2,x)
ans =
    (x*(2*x^2 + 9*x + 12))/6
>> expand(int(x^2 + 3*x + 2,x))
ans =
    x^3/3 + (3*x^2)/2 + 2*x

For definite integrals, we add the limits of integration to the arguments: int(s,var,a,b).

$\int_0^5 x^2 + 3\,x + 2 \, dx = \frac{535}{6}$

>> int(x^2 + 3*x + 2,x,0,5) % definite integral
ans =
    535/6

$\int_0^t x^2 + 3\,x + 2 \, dx = \frac{1}{3} t^3 + \frac{3}{2} t^2 + 2\,t$

>> expand(int(x^2 + 3*x + 2,x,0,t))
ans =
    t^3/3 + (3*t^2)/2 + 2*t

$\int_0^t e^{-x/2}\, dx = 2 - 2\,e^{-t/2}$

>> int(exp(-x/2),x,0,t)
ans =
    2 - 2*exp(-t/2)

$\int_0^\infty e^{-x/2}\, dx = 2$

>> int(exp(-x/2),x,0,Inf)  % definite to infinity
ans =
    2

When the symbolic math toolbox is not able to find an analytic solution, it may return a result that makes use of numerical integration methods. We also encountered the erf and erfc functions in Continuous Distributions.

>> p = int(exp(-x^2), -Inf, 2)
p =
(pi^(1/2)*(erf(2) + 1))/2

5.3.3. Symbolic Limits¶

The limit(s, h, k) function returns the limit of an expression, s as the variable h approaches k.

$\lim_{x \rightarrow 0} \frac{\sin(x)}{x} = 1$

$\lim_{x \rightarrow \infty} x\,e^{-x} = 0$

>> limit(sin(x)/x,x,0)
ans =
    1
>> limit(x*exp(-x),x,Inf)
ans =
    0

Another optional argument of 'left' or 'right' is useful when the expression has a discontinuity at $h = k$ . The expression $\frac{\abs{x - 3}}{x -3}$ has a discontinuity at $x = 3$ .

>> limit(abs(x-3)/(x-3),x,3,'right')
ans =
    1
>> limit(abs(x-3)/(x-3),x,3,'left')
ans =
    -1

The following computes $\frac{dy}{dx}$ for $y(x) = x^2$ using a limit. Recall the definition of a derivative:

$\frac{dy}{dx} = \lim_{h \rightarrow 0} \frac{y(x + h) - y(x)}{h}.$

>> syms h s x
>> s = ((x + h)^2 - x^2)/h;
>> limit(s,h,0)
ans =
    2*x

5.4. Symbolic Differential Equations¶

Although most students in this course have not yet studied the analytic math techniques for solving differential equations, we use the solutions to differential equations in all areas of physics and engineering. Differential equations are just expressions where one or more terms of the expression is the derivative of another variable. Any system that moves or changes in either time or space is defined by a differential equation.

The dsolve function returns the solution to differential equations. Note that previous versions of MATLAB used the special variables Dy and D2y for the first and second derivatives of $y$ . That notation is now deprecated. Instead use the diff function to define derivatives in the expression. As an example, consider the simple differential equation in (5.3), which is solved by the following MATLAB code.

(5.3)¶ $\begin{array}{l} \frac{dy(t)}{dt} = -a\,y(t) \\ \\ y(t) = C\,e^{-a\,t} \end{array}$

>> syms t y(t) a
>> y(t) = diff(y(t)) == -a*y(t);
>> Y(t) = dsolve(y(t))
Y(t) =
    C4*exp(-a*t)

The example above shows two important points about differential equations. First, the solution contains an exponential function of the irrational number $e$ . This is because $e^{a\,t}$ is the only function for which the derivative is a multiple of the original function. Specifically,

(5.4)¶ $\frac{d\,e^{a\,t}}{dt} = a\,e^{a\,t}.$

Thus, exponential functions of $e$ show up as the solution to any differential equation of growth or decay where there is a linear relationship between the rate of change and the function’s value.

Secondly, the solution to equation (5.4) contains a constant, C4. This is normal when an initial condition or boundary condition was not specified. With an initial condition, we can either solve for C4 algebraically, or use the condition as an input to dsolve. For example, let’s say that $Y(0) = 5$ . From our solution, $Y(0) = C4$ , thus $C4 = 5$ . The following code uses the SMT to solve a first order differential equation like (5.4) with exponential decay. The solution is plotted in Fig. 5.2.

Initial or boundary conditions may also be passed to dsolve as an extra argument. When there is more than one condition, as is common for second or higher order differential equations, the conditions are passed as an array. See the Throwing a Ball Up example for a second order differential equation example with two initial conditions.

>> syms y(t) a
>> eqn = diff(y(t)) == -a*y(t);
>> cond = y(0) == 5;
>> Y(t) = dsolve(eqn,cond)
Y(t) =
    5*exp(-a*t)

Note

As with function handles, the results computed by the SMT that are functions may be plotted with the fplot function.

>> a = 0.005;
>> y(t) = subs(Y(t))
y(t) =
    5*exp(-t/200)
>> fplot(y(t), [0 1000])

Fig. 5.2 The exponential decay solution to a differential equation of the same form as equation (5.4).