5. Using the Symbolic Math Toolbox

Before turning to algorithms for computing numeric solutions to problems, we need to point out that MATLAB can help find analytic solutions when such are available. Analytic solutions are always preferred over numeric solutions. When a computer is used to solve a math 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 algebraic equations describing the force of gravity on a moving ball are readily available on the web or in physics and dynamics textbooks. Let’s use the SMT to derive the equations and apply them.

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

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.

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

The equation above is a second order differential equation, which is our starting point. Things that move are described by differential equations. 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 with the syms statement. 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.

>> tm = solve(v(t) == 0)
tm =
    -V0/g
>> y_max = y(tm)
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.

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

We can give the variables numeric values and substitute (subs) the values into the equations. The vpa function gives decimal answers.

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

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

>> subs(tf(1))           % time until the ball hits the ground
ans =
    (5^(1/2)*598^(1/2))/49 + 50/49
>> vpa(subs(tf(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

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. Just be sure 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

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. In some cases, it may be desired to declare symbolic functions, such as g(x) below.

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

>> 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). For definite integrals, we add the limits of integration to the arguments: int(s,var,a,b).

\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

\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

5.3.3. Symbolic Limits

The int(s, h, k) function returns the limit of an expression, s as the variable h approaches k. Another optional argument of 'left' or 'right' is useful when the expression has a discontinuity at h = 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

The expression \frac{|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 class 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 to mean the first and second derivatives of y. That notation is now deprecated. Instead use the diff function to define any derivatives in the expression.

\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,

\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, decay, or resistance where there is a linear relationship between the derivative and the original function.

Secondly, the solution above 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.

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

The results computed by the SMT that are in the form of functions may be plotted with the fplot function. See Fzero for more on fplot.

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