Symbolic Toolbox and Related Functions
The Symbolic Math Toolkit is a Mathworks package that augments Matlab's existing functionality with the core Maple® symbolic kernel. With this package, you can solve and simplify systems of symbolic equations, find symbolic expressions for the inverse of a function, integrate, differentiate, take limits, and perform Taylor expansions, sums, variable precision arithmetic, or linear algebraic operations. If you have the extended symbolic toolkit, (which we do not discuss here), all of Maple's non-graphics packages are available.
Contents
There is more to this package than we can describe here. For more information, see the online documentation.
Working with Symbolic Variables
Symbolic variables are treated differently than regular variables in Matlab and must be created using the sym() or syms() functions.
syms A B lambda X Y Z a
Constant symbols can be defined too, which are not evaluated numerically.
r = sym(sqrt(2)/2) rr = r^2 t = sym(2/3) % Need to use sym here, not syms v = r + t w = r*2-3/t % notice that r is not evaluated q = sym(22/14 + 18/402) % add two fractions exactly
r = 2^(1/2)/2 rr = 1/2 t = 2/3 v = 2^(1/2)/2 + 2/3 w = 2^(1/2) - 9/2 q = 758/469
To convert a constant symbolic expression to a regular Matlab double value, use the double() function.
double(r) % numerically evaluate r
ans =
0.7071
We can build up more complicated symbolic expressions by defining functions of these variables. Note that these too are symbolic expressions, not function handles.
f = cos(X)
g = exp(X^2 - 2*X)
h = compose(f,g) % functional composition: f(g(X))
f = cos(X) g = exp(X^2 - 2*X) h = cos(exp(X^2 - 2*X))
Functions of multiple variables can also be created.
f = X^2 + Y^2 - 2*cos(X) g = sqrt(X^2 + Y^2) + Z^2 - lambda
f = X^2 - 2*cos(X) + Y^2 g = (X^2 + Y^2)^(1/2) - lambda + Z^2
The pretty() function tries to display a symbolic expression in a prettier way. It takes a bit of getting used to. Exponents, for example are printed on the line above, trying to mimic how you might write them by hand
pretty(g)
2 2 1/2 2
(X + Y ) - lambda + Z
You can convert an expression to latex as follows
latex(g)
ans =
\sqrt{X^2 + Y^2} - \mathrm{lambda} + Z^2
You can convert an expression to C code as follows
ccode(g)
ans = t0 = -lambda+sqrt(X*X+Y*Y)+Z*Z;
Symbolic matrices are created in much the same way numeric matrices are.
syms E F G H I J mat = [E F G ; H I J]
mat = [ E, F, G] [ H, I, J]
We will make use of these in the sections to come.
The subs() function can be used to substitute one value for another, including a numeric value for symbolic one.
subs(f,X,3) % substitute 3 for X in f
ans = Y^2 - 2*cos(3) + 9
Once all of the symbolic variables are numeric, the result is numerically evaluated. We can prevent this by substituting sym(3) and sym(10) instead.
subs(f,{X,Y},{3,10}) % substitute 3 for X, 10 for Y in f
ans = 110.9800
subs(f,X,lambda-Y) % substitute (lambda - Y) for X in f
ans = Y^2 - 2*cos(lambda - Y) + (Y - lambda)^2
Many operations on symbolic expressions are ambiguous unless the independent variable is specified. If this is not given explicitly, Matlab chooses the variable closest in alphabetical order to x, (ties broken in favor of the end of the alphabet). You can see how Matlab will order the variables in an expression with the findsym() function.
findsym(f)
ans = X,Y
Basic Algebra
There are a number of functions that we can use to perform basic high school algebra that might be tedious or error prone to do by hand. We have already seen compose(), which belongs on this list.
S = (((2*X + 8*X^2)/(2*X)+(2*X^2 - 2*X + X)*(X+2)-(2*X^2 + 2)/(4*X^2 - 2*X^2)*(X+2))/(X+1))+1; pretty(S);
/ 2
| 2 8 X + 2 X
1 - | (X + 2) (X - 2 X ) - ---------- +
| 2 X
\
2 \
(2 X + 2) (X + 2) |
------------------ | / (X + 1)
2 |
2 X /
We can then factor or expand S.
Sf = factor(S); % factor S
pretty(Sf);
5 4 3
2 X + 3 X + 2 X - X - 2
--------------------------
2
X (X + 1)
Se = expand(S); % expand s
pretty(Se);
2 3
3 X 1 2 X 2 X
----- - ----- + ----- - ------- + ----- -
X + 1 X + 1 X + 1 3 2 X + 1
X + X
1
------ + 1
2
X + X
The simple() and simplify() functions can be used to try and find the simplest representation: simple tends to do better with trigonometric expressions. Surprisingly, we can sometimes get simpler expressions by applying the function multiple times.
T = 1/2*(2*tan(1/2*X)/(1+tan(1/2*X)^2) *... (1-tan(1/2*Y)^2)/(1+tan(1/2*Y)^2)-2 *... (1-tan(1/2*X)^2)/(1+tan(1/2*X)^2) *... tan(1/2*Y)/(1+tan(1/2*Y)^2))/tan(1/2*X-1/2*Y) *... (1+tan(1/2*X-1/2*Y)^2); Tsimple = simple(T) TverySimple = simple(simple(T))
Tsimple = 1 TverySimple = 1
If you call simple() without saving the output, i.e. just simple(T) as opposed to t = simple(T), it displays many equivalent expressions.
The collect() function can be used to collect like terms. Here we collect all of the Y variables together and find we have (3X - 3Z) of them
t = (X + Y)*(2*X-3*Z)+Z+X*Y-(3*X+4); pretty(t);
Z - 3 X + X Y + (X + Y) (2 X - 3 Z) - 4
t = collect(t,Y); pretty(t);
(3 X - 3 Z) Y + Z - 3 X + X (2 X - 3 Z) - 4
Function Inverse
We can find the inverse of a function, (if it exists) with finverse().
f = 2*sin(cos(3*log(X)/4))+1 finv = finverse(f)
f = 2*sin(cos((3*log(X))/4)) + 1 Warning: finverse(2*sin(cos((3*ln(X))/4)) + 1) cannot be found. finv = [ empty sym ]
Solving Symbolic Equations
The solve() function can be used to solve systems of equations, symbolically. Its inputs are either symbolic expressions or strings, with each equation separated by a comma followed by the variable or variables you wish to solve for. The output is a struct with a field for each variable.
S = solve('k = a/(a-b+p)','p = k*(1+k)^2/(2*a + b+ 1)','a','b')
S =
a: [1x1 sym]
b: [1x1 sym]
pretty(S.a)
4 3 2 2
k + 2 k + k - k p - k p
- ---------------------------
p - 3 k p
pretty(S.b)
4 3 2 2
k + k - k + 2 k p - k p - k + p
- -----------------------------------
p - 3 k p
If the equations do not contain an equals sign, they are assumed equal to 0.
pretty(solve('a*X^2+b*X+c','X'))
+- -+ | 2 1/2 | | b + (b - 4 a c) | | - ------------------- | | 2 a | | | | 2 1/2 | | b - (b - 4 a c) | | - ------------------- | | 2 a | +- -+
S = solve('X = 2*Y-1','Y=3*Z','Z=X+Y','X','Y','Z'); S.X S.Y S.Z
ans = -1/4 ans = 3/8 ans = 1/8
Limits
We can take limits of functions, (two-sided, as well as left and right) as they approach specific values, (symbolic or numeric) or as they tend to inf or -inf.
limit((X-1)/(sqrt(X)-1),1) % limit as X --> 2 limit(X^X/(exp(X)^log(X)),inf) % limit as X --> inf limit(cos(X)-2*3^X,X,a) % limit as X --> a limit(cos(X*Y)^X-2*3^Y,Y,0) % limit as Y --> 0 limit(exp(Y-X),Y,-inf) % limit as Y --> -inf
ans = 2 ans = 1 ans = cos(a) - 2*3^a ans = -1 ans = 0
f = 1/(1 + 2^(-1/X)); limit(f,0) % two sided limit as X --> 0 (D.N.E.) limit(f,X,0,'left') % left sided limit: X-->0- limit(f,X,0,'right') % right sided limit: X-->0+
ans = NaN ans = 0 ans = 1
Symbolic Sums
The symsum() function performs symbolic sums. We specify the summand w.r.t. the indexing variable, followed by the start and end indices. The end index can be inf, in which case an infinite sum is performed.
symsum(2^(-X),0,inf) symsum((-1)^(X+1)*(1/X),1,inf) symsum(1/(X^2),1,inf) symsum((-1)^(X+1)*(1/(2*X-1)),1,inf) symsum(1/X,1,10)
ans = 2 ans = log(2) ans = pi^2/6 ans = pi/4 ans = 7381/2520
Taylor Series Expansions
We can perform a Taylor expansion of a function. We specify the function, and optionally, the degree of the expansion, (default is 6), and a symbol or numeric value about which the expansion is performed.
taylor(sin(X)) % 6th order Taylor expansion of sin(X) taylor(exp(X)) % 6th order Taylor expansion of exp(X)
ans = X^5/120 - X^3/6 + X ans = X^5/120 + X^4/24 + X^3/6 + X^2/2 + X + 1
syms A pretty(taylor(sin(X),10,A)) % 10th order expansion about the point A
3
cos(A) (A - X)
sin(A) - cos(A) (A - X) + --------------- -
6
5 7
cos(A) (A - X) cos(A) (A - X)
--------------- + --------------- -
120 5040
9 2
cos(A) (A - X) sin(A) (A - X)
--------------- - --------------- +
362880 2
4 6
sin(A) (A - X) sin(A) (A - X)
--------------- - --------------- +
24 720
8
sin(A) (A - X)
---------------
40320
pretty(taylor(int(normpdf(X)),10)); % 10th order expansion of the integral of the normal distribution
9 7
17592186044416 X 140737488355328 X
------------------ - ------------------ +
152399477208391047 118532926717637481
5 3
281474976710656 X 1125899906842624 X
------------------ - ------------------- +
28222125408961305 16933275245376783
2251799813685248 X
------------------
5644425081792261
Differentiation
The diff() function performs differentiation.
diff(log(X)) % differentiate with respect to X diff(log(X),3) % take the 3rd derivative of log(X) w.r.t. X diff(Y*sin(X)+2*X^Y,'Y') % differentiate w.r.t Y
ans = 1/X ans = 2/X^3 ans = sin(X) + 2*X^Y*log(X)
We can combine diff() with solve() to find maxima and minima
solve((diff(sqrt(log(3*X^3-2*X)))))
ans = -2^(1/2)/3 2^(1/2)/3
Here we attempt to find the MLE for the product of N univariate Gaussians. Unfortunately, we have to specify a value for N and proceed via induction.
syms mu s2 pi normconst = sqrt(2*pi*s2); % s2 for sigma^2 NormPDF = (1/normconst)*exp((-(X-mu)^2)/(2*s2)); % single Gaussian distribution prodNorm = 1; N = 5; % we will start with N=5 for i=1:N syms(['X',num2str(i)]); prodNorm = prodNorm*subs(NormPDF,X,['X',num2str(i)]); % take product of N Gaussians end muMLE = solve(diff(prodNorm,mu),mu); % MLE w.r.t. mu s2MLE = solve(diff(prodNorm,s2),s2); % MLE w.r.t. sigma^2 pretty(muMLE)
X1 X2 X3 X4 X5 -- + -- + -- + -- + -- 5 5 5 5 5
We can easily spot the pattern: (1/N)sum(X)
Unfortunately, s2MLE is not quite as pretty.
pretty(s2MLE)
2 2 2
X1 2 X1 mu X2 2 X2 mu X3 2 X3 mu
--- - ------- + --- - ------- + --- - ------- +
5 5 5 5 5 5
2 2
X4 2 X4 mu X5 2 X5 mu 2
--- - ------- + --- - ------- + mu
5 5 5 5
Integration
The int() function can perform definite and indefinite integration.
int(log(X)) %indefinite integral w.r.t X int(sin(X)^2) %indefinite integral w.r.t X int(int(log(X+Y),'X'),'Y') %indefinite integral w.r.t Y
ans = X*(log(X) - 1) ans = X/2 - sin(2*X)/4 ans = ((X + Y)^2*(2*log(X + Y) - 3))/4
result = int(2*X*log(X),1,10) % definite integral double(result) % convert exact to a double
result = 100*log(10) - 99/2 ans = 180.7585
Multivariate Calculus
Here we take the Jacobian of f with respect to v, where both are vectors.
syms r theta v = [r theta]; x1 = r*cos(theta); x2 = r*sin(theta); f = [x1,x2]; J = jacobian(f,v)
J = [ cos(theta), -r*sin(theta)] [ sin(theta), r*cos(theta)]
We can calculate the determinant exactly.
detJ = (det(J))
detJ = r*cos(theta)^2 + r*sin(theta)^2
Lets simplify the above expression.
detJ = simplify(detJ)
detJ = r
We can also use this function to compute a gradient, so long as our function f is scalar valued.
syms X Y Z f = 2*sin(X) + cos(Y)^2 -3*log(Z) J = jacobian(f,[X,Y,Z]); pretty(J);
f = cos(Y)^2 - 3*log(Z) + 2*sin(X) +- -+ | 3 | | 2 cos(X), (-2) cos(Y) sin(Y), - - | | Z | +- -+
Linear Algebra
We can create symbolic matrices as we saw above, and perform various operations. Constant symbols can be used to perform computations exactly.
clear syms a b c d e f g h i j k l m n o p A = [a b ; c d ; e f]; B = [g h i ; j k l]; C = A*B
C = [ a*g + b*j, a*h + b*k, a*i + b*l] [ c*g + d*j, c*h + d*k, c*i + d*l] [ e*g + f*j, e*h + f*k, e*i + f*l]
D = [a b ; c d]; inv(D)
ans = [ d/(a*d - b*c), -b/(a*d - b*c)] [ -c/(a*d - b*c), a/(a*d - b*c)]
det(D)
ans = a*d - b*c
eig(D)
ans = a/2 + d/2 - (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2 a/2 + d/2 + (a^2 - 2*a*d + d^2 + 4*b*c)^(1/2)/2
E = sym([1/3 4/15 2/5 ; 2/3 1/3 0; 1/2 1/4 1/4]) E2 = E^2
E = [ 1/3, 4/15, 2/5] [ 2/3, 1/3, 0] [ 1/2, 1/4, 1/4] E2 = [ 22/45, 5/18, 7/30] [ 4/9, 13/45, 4/15] [ 11/24, 67/240, 21/80]
det(E)
ans = -1/60
inv(E)
ans = [ -5, -2, 8] [ 10, 7, -16] [ 0, -3, 4]
Here we try and compute the SVD exactly but to no avail. We can find a rational approximation to the entries, however, by using the rats() function, which we discuss in more depth below.
[U S V] = svd(E) % compute the SVD of E
U = rats(double(U))
S = rats(double(S))
V = rats(double(V))
U =
[ 0.49109823962669798664999567078201, -0.72406791233697925247981479395538, -0.48430174205708454355751477479149]
[ 0.65982853521141748292248273038817, 0.67217592578665185464712183886594, -0.33586579003169211900746006123949]
[ 0.56872561324096659876202200872008, -0.15461301082646966779516726932017, 0.80786508386416877686874947261274]
S =
[ 1.0706292924174921593709666623039, 0, 0]
[ 0, 0.35296678714206734207932738121935, 0]
[ 0, 0, 0.044103777275924655598566463054976]
V =
[ 0.82937008753546424369310341285868, 0.36676200815867437795437640947586, 0.4214627946483922404785452801969]
[ 0.46055509823988214272694229122357, -0.021756685559871156482933444224287, -0.88736443929126515473068325503709]
[ 0.31628193022471088744597118209586, -0.93006036148460029020028489223978, 0.18695845690544664601643216465309]
U =
331/674 -349/482 -108/223
225/341 244/363 -221/658
571/1004 -62/401 185/229
S =
1531/1430 0 0
0 809/2292 0
0 0 46/1043
V =
593/715 128/349 161/382
216/469 -27/1241 -323/364
68/215 -625/672 43/230
Symbolic and non-symbolic expressions can be combined.
[a b ; c d]*[2 0 ; 3 5]
ans = [ 2*a + 3*b, 5*b] [ 2*c + 3*d, 5*d]
Here we test the correctness of the Matrix Inversion Lemma. We define A,B,C, and D to be arbitrary square matrices. Unfortunately, their size is fixed and so the test is not perfectly general.
Keep in mind that two symbolic equations are only considered equal if Matlab is able to find identical representations of each. Its searches are not exhaustive and as such, it might claim that two expressions are unequal, when in fact they are, but not conversely.
syms a b c d e f g h i j k l m n o p A = [ a b ; c d]; B= [e f ; g h]; C = [ i j ; k l]; D = [m n ; o p]; lhs = inv([A B ; C D]); rhs11 = inv(A) + inv(A)*B*inv(D-C*inv(A)*B)*C*inv(A); rhs12 = -inv(A)*B*inv(D - C*inv(A)*B); rhs21 = -inv(D-C*inv(A)*B)*C*inv(A); rhs22 = inv(D-C*inv(A)*B); rhs = [rhs11 rhs12 ; rhs21 rhs22]; test = simplify(lhs) == simplify(rhs)
test =
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
Variable Precision Arithmetic
The vpa() function lets us perform variable precision arithmetic. The default is 32 digits, which can be changed for the current session with the digits() function.
PI = vpa(pi,80) % pi to 80 digits s = vpa(sin(sqrt(2)/2),40) % calculate to 40 digits
PI = 3.141592653589793115997963468544185161590576171875 s = 0.6496369390800624810111685292213223874569
Rational Fraction Approximation
As the title suggests, we an approximate any real number as a rational fraction, (whose numerator and denominator are relatively prime). This is actually a Matlab function, (not from the symbolic toolbox).
q = rats((22/14) + (18/402)) e = rats(exp(1)) simplified = rats(940/40) val = rats(9.352422118484)
q =
758/469
e =
1457/536
simplified =
47/2
val =
2123/227
clear