# 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 to , but would leave numbers like , , and 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 . The velocity of the ball is the derivative of the position,

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.

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 ; and the initial velocity from the throw is .

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 is 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 my 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)`

.

```
>> 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(x^2 + 3*x + 2,x,0,5) % definite integral
ans =
535/6
```

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

```
>> int(exp(-x/2),x,0,t)
ans =
2 - 2*exp(-t/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
.

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

The expression has a discontinuity at .

```
>> 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 for using a limit. Recall the definition of a derivative:

```
>> 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 . That notation
is now deprecated. Instead use the `diff`

function to define any
derivatives in the expression.

```
>> 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 . This is because is the only function for which the derivative is a multiple of the original function. Specifically,

Thus, exponential functions of 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 . From our solution,
, thus .

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 Zero Finding for more on
`fplot`

.

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