5. Using the Symbolic Math Toolbox

We now take a break from numerical algorithms and focus on analytic solutions to mathematical problems. We seek analytical solutions in mathematics courses and 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 mathematical 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.

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 figure 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. 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.

A diagram showing the path of a ball thrown upwards

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\).

Let \(y(t)\) represent the vertical position of the ball with respect to time. 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 negative constant due to gravity. 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/sec2, or \(g = -32\) feet/sec2

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 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  %  second order derivative
diffeq =
    diff(y(t), t, 2) == g
>> v(t) = diff(y(t),t)          %  first order derivative
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. We want the first answer (t_final(1) because it is the positive number.

>> 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 the ground, 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 forms 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.

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)

>> %  intersection of
>> %  two functions
>> 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)
>> %  non-polynomial
>> y(t) = cos(pi*t)
y(t) =
    cos(pi*t)
>> solve(y(t), t)
ans =
    1/2

>> y(t) = 1 + t^2;
>> %  finds complex roots
>> solve(y(t), t)
ans =
   -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 decimal approximations. The vpa function converts symbolic expressions to decimal values. 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
>> %  2nd derivative
>> diff(x^2 + 3*x + 2, x, 2)
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 \left(x^2 + 3\,x + 2\right) \, 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 \left(x^2 + 3\,x + 2\right) \, dx = \frac{535}{6}\]
>> int(x^2 + 3*x + 2, x, 0, 5) % definite integral
ans =
    535/6

The limits of integration may be another variable.

\[\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)

Infinity may also be an integration limit.

\[\int_0^\infty e^{-x/2}\, dx = 2\]
>> int(exp(-x/2), x, 0, Inf)  % zero to infinity
ans =
    2

When the SMT cannot find an analytic solution, it may return a result that uses 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

We use the solutions to differential equations in all areas of physics and engineering. Differential equations are expressions where one or more terms of the expression are the derivative of a variable. A differential equation defines any system that moves or changes in either time or space.

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. For example, consider the simple differential equation in equation (5.4), which is solved by the following MATLAB code.

(5.4)\[\begin{split}\begin{array}{l} \frac{dy(t)}{dt} = -a\,y(t) \\ \\ y(t) = C\,e^{-a\,t} \end{array}\end{split}\]
>> 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\) because \(e^{a\,t}\) is a function for which the derivative is a multiple of the original function. Specifically,

(5.5)\[\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.5) contains a constant, C4, which 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 equation (5.4) with exponential decay. The solution is plotted in figure Fig. 5.2.

Initial or boundary conditions may also be passed to dsolve. 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)

>> a = 0.005;
>> y(t) = subs(Y(t))
y(t) =
    5*exp(-t/200)
>> fplot(y(t), [0 1000])
The exponential decay solution to a differential equation

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