Symbolic Toolbox and Related Functions in Matlab
by Matt Dunham
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.
Useful Functions
sym, syms, double, pretty, subs, findsym, compose, factor, expand, collect, subexpr, finverse, simplify, simple, solve, limit, symsum, taylor, diff, int, jacobian, inv, det, eig, svd, vpa, rats,
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 = sqrt(1/2) rr = 1/2 t = 2/3 v = 1/2*2^(1/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+Y^2-2*cos(X) g = (X^2+Y^2)^(1/2)+Z^2-lambda
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, (or how latex would typeset them).
pretty(g)
2 2 1/2 2 (X + Y ) + Z - lambda
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 = 9+Y^2-2*cos(3)
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 = (lambda-Y)^2+Y^2-2*cos(lambda-Y)
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, YBasic 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 2 X + 8 X 2 (2 X + 2) (X + 2) 1/2 ---------- + (2 X - X) (X + 2) - 1/2 ------------------ X 2 X ------------------------------------------------------------ + 1 X + 1
We can then factor or expand S.
Sf = factor(S); % factor S pretty(Sf);
3 5 4 2 X + 2 X + 3 X - X - 2 -------------------------- 2 X (X + 1)
Se = expand(S); % expand s pretty(Se);
3 2 1 X X X 1 2 - ----- + ----- + 2 ----- + 3 ----- - --------- - ---------- + 1 X + 1 X + 1 X + 1 X + 1 (X + 1) X 2 (X + 1) 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/2*(sin(Y)*cos(X)-sin(X)*cos(Y))/cos(1/2*X-1/2*Y)/sin(1/2*X-1/2*Y) 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);
(X + Y) (2 X - 3 Z) + Z + X Y - 3 X - 4
t = collect(t,Y); pretty(t);
(3 X - 3 Z) Y + X (2 X - 3 Z) - 4 + Z - 3 X
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/4*log(X)))+1 finv = exp(4/3*acos(asin(-1/2+1/2*X)))
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)
3 2 2 k (k - p + k + 2 k - p ) -------------------------- p (3 k - 1)
pretty(S.b)
4 2 3 2 k - k p - k + k + 2 k p + p - k ----------------------------------- p (3 k - 1)
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) ] [- 1/2 -------------------] [ a ] [ ] [ 2 1/2] [ b + (b - 4 a c) ] [- 1/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 = 1/6*pi^2 ans = 1/4*pi 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-1/6*X^3+1/120*X^5 ans = 1+X+1/2*X^2+1/6*X^3+1/24*X^4+1/120*X^5
syms A pretty(taylor(sin(X),10,A)) % 10th order expansion about the point A
2 3 sin(A) + cos(A) (X - A) - 1/2 sin(A) (X - A) - 1/6 cos(A) (X - A) 4 5 + 1/24 sin(A) (X - A) + 1/120 cos(A) (X - A) 6 7 - 1/720 sin(A) (X - A) - 1/5040 cos(A) (X - A) 8 9 + 1/40320 sin(A) (X - A) + 1/362880 cos(A) (X - A)
pretty(taylor(int(normpdf(X)),10)); % 10th order expansion of the integral of the normal distribution
2251799813685248 1125899906842624 3 281474976710656 5 ---------------- X - ----------------- X + ----------------- X 5644425081792261 16933275245376783 28222125408961305 140737488355328 7 17592186044416 9 - ------------------ X + ------------------ X 118532926717637481 152399477208391047
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 = 1/3*2^(1/2) -1/3*2^(1/2)
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)
1/5 X1 + 1/5 X2 + 1/5 X3 + 1/5 X4 + 1/5 X5
We can easily spot the pattern: (1/N)sum(X)
Unfortunately, s2MLE is not quite as pretty.
pretty(s2MLE)
2 2 2 2 1/5 X1 - 2/5 X1 mu + mu + 1/5 X2 - 2/5 X2 mu + 1/5 X3 - 2/5 X3 mu 2 2 + 1/5 X4 - 2/5 X4 mu + 1/5 X5 - 2/5 X5 mu
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)-X ans = -1/2*sin(X)*cos(X)+1/2*X ans = 1/2*log(X+Y)*(X+Y)^2-1/4*X^2-3/2*X*Y-3/4*Y^2
result = int(2*X*log(X),1,10) % definite integral double(result) % convert exact to a double
result = -99/2+100*log(2)+100*log(5) 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 = cos(theta)^2*r+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 = 2*sin(X)+cos(Y)^2-3*log(Z) [2 cos(X) -2 cos(Y) sin(Y) - 3/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*ceig(D)
ans = 1/2*a+1/2*d+1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2) 1/2*a+1/2*d-1/2*(a^2-2*a*d+d^2+4*b*c)^(1/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 = [ -.49109823962669798664999567078201, .72406791233697925247981479395543, -.48430174205708454355751477479149] [ -.65982853521141748292248273038823, -.67217592578665185464712183886597, -.33586579003169211900746006123953] [ -.56872561324096659876202200872013, .15461301082646966779516726932019, .80786508386416877686874947261276] S = [ 1.0706292924174921593709666623038, 0, 0] [ 0, .35296678714206734207932738121935, 0] [ 0, 0, .44103777275924655598566463054979e-1] V = [ -.82937008753546424369310341285867, -.36676200815867437795437640947584, .42146279464839224047854528019686] [ -.46055509823988214272694229122354, .21756685559871156482933444224305e-1, -.88736443929126515473068325503707] [ -.31628193022471088744597118209584, .93006036148460029020028489223976, .18695845690544664601643216465313] 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 1Variable 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.1415926535897932384626433832795028841971693993751058209749445923078164062862090 s = .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/227clearPopularity: 1% [?]
Automatic Numerical Differentiation
Computing the differential properties of function can often be an heavy task. On Maltab file exchange there is an advanced tool that can make our life easier. This tool is is Automatic Numerical Differentiation from John D’Errico.
It allows for Numerical derivative of an analytically supplied function, also gradient, Jacobian & Hessian.
The author provides the following brief description:
The DERIVESTsuite now provides automatic numerical differentiation for both scalar and vector valued functions. Tools for derivatives (up to 4th order) of a scalar function are provided, as well as the gradient vector, directional derivative, Jacobian matrix, and Hessian matrix. Error estimates are provided for all tools.
DERIVEST provides automatic numerical differentiation (up to the fourth derivative) of a user supplied function, much as quad does for integration. It is semi-intelligent, trying to use that step size which minimizes its estimate of the uncertainty in the derivative.
High order methods are used, although full control is provided to the user when you want it. You can direct the order of the method to be used, the general class of difference method employed (forward, backward, or central differences), the number of terms employed in its generalized Romberg acceleration, step sizes, etc.
Although you can not provide a user supplied tolerance, DERIVEST does return an estimate of its uncertainty in the final result.
For example, the derivative of exp(x), at x=1 is exp(1)==2.71828182845905. DERIVEST comes pretty close.
[d,err]=derivest(@(x) exp(x),1) d = 2.71828182845904 err = 1.02015503167879e-14
Many more examples can be found at following links:
- Gradient of the Rosenbrock function at [1,1], the global minimizer
- The Hessian matrix at the minimizer should be positive definite
- Gradient estimation using gradest – a function of 5 variables
- Simple Hessian matrix of a problem with 3 independent variables
- A semi-definite Hessian matrix
- Directional derivative of the Rosenbrock function at the solution
- Directional derivative at other locations
- Jacobian matrix of a scalar function is just the gradient
- Jacobian matrix of a linear system will reduce to the design matrix
- The jacobian matrix of a nonlinear transformation of variables
- derivative of exp(x), at x == 0
- DERIVEST can also use an inline function
- Higher order derivatives (second derivative)
- Higher order derivatives (third derivative)
- Higher order derivatives (up to the fourth derivative)
- Evaluate the indicated (default = first) derivative at multiple points
- Specify the step size (default stepsize = 0.1)
- Provide other parameters via an anonymous function
- The second derivative should be positive at a minimizer.
- Compute the numerical gradient vector of a 2-d function
- Compute the numerical Laplacian function of a 2-d function
- Compute the derivative of a function using a central difference scheme
- Compute the derivative of a function using a forward difference scheme
- Compute the derivative of a function using a backward difference scheme
- Although a central rule may put some samples in the wrong places, it may still succeed
- But forcing the use of a one-sided rule may be smart anyway
- Control the behavior of DERIVEST – forward 2nd order method, with only 1 Romberg term
- Functions should be vectorized for speed, but its not always easy to do.
Popularity: 1% [?]
MatlabBGL: a graph library for Matlab
Filed under: Computational Geometry, Mathematics, Numerical Analysis
The most complete graph library for Matlab offered by David Gleich.
The MatlabBGL library fills a hole in Matlab’s suite of algorithms. Namely, it provides a rich set of algorithms to work with graphs, as in graph theory graphs. The MatlabBGL package uses Matlab’s native sparse matrix type as a graph and provides algorithms that work
The algorithms included are
- Searching: breadth first search,depth first search, and astar (A*) search
- Shortest Path Algorithms: Dijkstra’s algorithm, the Bellman-Ford algorithm, Johnson’s algorithm, and the Floyd-Warshall algorithm.
- Minimum Spanning Trees: Prim’s algorithm and Kruskal’s algorithm.
- Components: strongly connected components and biconnected components (and articulation points).
- Flow Algorithms: Goldberg’s push-relabel maximum-flow minimum-cut algorithm.
- Statistics: Betweenness Centrality, Clustering Coefficients, and Edge Centrality
- Graph Creation: Erdos Reyni (Gnp) Graph, Cycle Graph, Wheel Graph, Star Graph
- Planar Graphs: Boyer-Myrvold planarity testing, Chrobak-Payne straight line drawing
- Graph Layout: force directed layout, spring based layout, topology filling layout
Additional informations can be found on the MatlabBGL official website, here is a directory:
MatlabBGL
Popularity: 1% [?]
Translator
Follow Us
Video Posts
