Store 
| Products & Services | Solutions | Academia | Support | User Community | Company |
 |
|
|
|
|
|
 |
Download Product Updates | | |  |
Get Pricing | | | |
Trial Software |
| R2010b Documentation → Symbolic Math Toolbox |
| Products & Services | Solutions | Academia | Support | User Community | Company |
|
|
|
|
|
|
|
Download Product Updates | | | |
Get Pricing | | | |
Trial Software |
| R2010b Documentation → Symbolic Math Toolbox |
| Contents | Index |
| Learn more about Symbolic Math Toolbox |
| On this page… |
|---|
Basic algebraic operations on symbolic objects are the same as operations on MATLAB objects of class double. This is illustrated in the following example.
The Givens transformation produces a plane rotation through the angle t. The statements
syms t; G = [cos(t) sin(t); -sin(t) cos(t)]
create this transformation matrix.
G = [ cos(t), sin(t)] [ -sin(t), cos(t)]
Applying the Givens transformation twice should simply be a rotation through twice the angle. The corresponding matrix can be computed by multiplying G by itself or by raising G to the second power. Both
A = G*G
and
A = G^2
produce
A = [ cos(t)^2 - sin(t)^2, 2*cos(t)*sin(t)] [ -2*cos(t)*sin(t), cos(t)^2 - sin(t)^2]
The simple function
A = simple(A)
uses a trigonometric identity to return the expected form by trying several different identities and picking the one that produces the shortest representation.
A = [ cos(2*t), sin(2*t)] [ -sin(2*t), cos(2*t)]
The Givens rotation is an orthogonal matrix, so its transpose is its inverse. Confirming this by
I = G.' *G
which produces
I = [ cos(t)^2 + sin(t)^2, 0] [ 0, cos(t)^2 + sin(t)^2]
and then
I = simple(I)
I = [ 1, 0] [ 0, 1]
The following examples show how to do several basic linear algebraic operations using Symbolic Math Toolbox software.
The command
H = hilb(3)
generates the 3-by-3 Hilbert matrix. With format short, MATLAB prints
H =
1.0000 0.5000 0.3333
0.5000 0.3333 0.2500
0.3333 0.2500 0.2000The computed elements of H are floating-point numbers that are the ratios of small integers. Indeed, H is a MATLAB array of class double. Converting H to a symbolic matrix
H = sym(H)
gives
H = [ 1, 1/2, 1/3] [ 1/2, 1/3, 1/4] [ 1/3, 1/4, 1/5]
This allows subsequent symbolic operations on H to produce results that correspond to the infinitely precise Hilbert matrix, sym(hilb(3)), not its floating-point approximation, hilb(3). Therefore,
inv(H)
produces
ans = [ 9, -36, 30] [ -36, 192, -180] [ 30, -180, 180]
and
det(H)
yields
ans = 1/2160
You can use the backslash operator to solve a system of simultaneous linear equations. For example, the commands
% Solve Hx = b b = [1; 1; 1]; x = H\b
produce the solution
x = 3 -24 30
All three of these results, the inverse, the determinant, and the solution to the linear system, are the exact results corresponding to the infinitely precise, rational, Hilbert matrix. On the other hand, using digits(16), the command
digits(16); V = vpa(hilb(3))
returns
V = [ 1.0, 0.5, 0.3333333333333333] [ 0.5, 0.3333333333333333, 0.25] [ 0.3333333333333333, 0.25, 0.2000000000007276]
The decimal points in the representation of the individual elements are the signal to use variable-precision arithmetic. The result of each arithmetic operation is rounded to 16 significant decimal digits. When inverting the matrix, these errors are magnified by the matrix condition number, which for hilb(3) is about 500. Consequently,
inv(V)
which returns
ans = [ 9.000000000000061, -36.00000000000032, 30.0000000000003] [ -36.00000000000032, 192.0000000000017, -180.0000000000015] [ 30.0000000000003, -180.0000000000015, 180.0000000000014]
shows the loss of two digits. So does
1/det(V)
which gives
ans = 2160.000000000018
and
V\b
which is
ans = 3.000000000000041 -24.00000000000021 30.00000000000019
Since H is nonsingular, calculating the null space of H with the command
null(H)
returns an empty matrix:
ans = [ empty sym ]
Calculating the column space of H with
colspace(H)
returns a permutation of the identity matrix:
ans = [ 1, 0, 0] [ 0, 1, 0] [ 0, 0, 1]
A more interesting example, which the following code shows, is to find a value s for H(1,1) that makes H singular. The commands
syms s H(1,1) = s Z = det(H) sol = solve(Z)
produce
H = [ s, 1/2, 1/3] [ 1/2, 1/3, 1/4] [ 1/3, 1/4, 1/5] Z = s/240 - 1/270 sol = 8/9
Then
H = subs(H, s, sol)
substitutes the computed value of sol for s in H to give
H = [ 8/9, 1/2, 1/3] [ 1/2, 1/3, 1/4] [ 1/3, 1/4, 1/5]
Now, the command
det(H)
returns
ans = 0
and
inv(H)
produces the message
ans = FAIL
because H is singular. For this matrix, null space and column space are nontrivial:
Z = null(H) C = colspace(H)
Z =
3/10
-6/5
1
C =
[ 1, 0]
[ 0, 1]
[ -3/10, 6/5]It should be pointed out that even though H is singular, vpa(H) is not. For any integer value d, setting digits(d), and then computing inv(vpa(H)) results in an inverse with elements on the order of 10^d.
The symbolic eigenvalues of a square matrix A or the symbolic eigenvalues and eigenvectors of A are computed, respectively, using the commands E = eig(A) and [V,E] = eig(A).
The variable-precision counterparts areE = eig(vpa(A)) and [V,E] = eig(vpa(A)).
The eigenvalues of A are the zeros of the characteristic polynomial of A, det(A-x*I), which is computed by poly(A).
The matrix H from the last section provides the first example:
H = sym([8/9 1/2 1/3; 1/2 1/3 1/4; 1/3 1/4 1/5])
H = [ 8/9, 1/2, 1/3] [ 1/2, 1/3, 1/4] [ 1/3, 1/4, 1/5]
The matrix is singular, so one of its eigenvalues must be zero. The statement
[T,E] = eig(H)
produces the matrices T and E. The columns of T are the eigenvectors of H and the diagonal elements of E are the eigenvalues of H:
T = [ 218/285 - (4*12589^(1/2))/285, (4*12589^(1/2))/285 + 218/285, 3/10] [ 292/285 - 12589^(1/2)/285, 12589^(1/2)/285 + 292/285, -6/5] [ 1, 1, 1] E = [ 32/45 - 12589^(1/2)/180, 0, 0] [ 0, 12589^(1/2)/180 + 32/45, 0] [ 0, 0, 0]
It may be easier to understand the structure of the matrices of eigenvectors, T, and eigenvalues, E, if you convert T and E to decimal notation. To do so, proceed as follows. The commands
Td = double(T) Ed = double(E)
return
Td =
-0.8098 2.3397 0.3000
0.6309 1.4182 -1.2000
1.0000 1.0000 1.0000
Ed =
0.0878 0 0
0 1.3344 0
0 0 0The first eigenvalue is zero. The corresponding eigenvector
(the first column of Td) is the same as the basis
for the null space found in the last section. The other two eigenvalues
are the result of applying the quadratic formula to
which is the quadratic factor
in factor(poly(H)):
syms x g = simple(factor(poly(H))/x); solve(g)
ans = 32/45 - 12589^(1/2)/180 12589^(1/2)/180 + 32/45
Closed form symbolic expressions for the eigenvalues are possible only when the characteristic polynomial can be expressed as a product of rational polynomials of degree four or less. The Rosser matrix is a classic numerical analysis test matrix that illustrates this requirement. The statement
R = sym(rosser)
generates
R = [ 611, 196, -192, 407, -8, -52, -49, 29] [ 196, 899, 113, -192, -71, -43, -8, -44] [ -192, 113, 899, 196, 61, 49, 8, 52] [ 407, -192, 196, 611, 8, 44, 59, -23] [ -8, -71, 61, 8, 411, -599, 208, 208] [ -52, -43, 49, 44, -599, 411, 208, 208] [ -49, -8, 8, 59, 208, 208, 99, -911] [ 29, -44, 52, -23, 208, 208, -911, 99]
The commands
p = poly(R); pretty(factor(p))
produce
2 2 2
x (x - 1020) (x - 1040500) (x - 1020 x + 100) (x - 1000)The characteristic polynomial (of degree 8) factors nicely into the product of two linear terms and three quadratic terms. You can see immediately that four of the eigenvalues are 0, 1020, and a double root at 1000. The other four roots are obtained from the remaining quadratics. Use
eig(R)
to find all these values
ans =
0
1000
1000
1020
510 - 100*26^(1/2)
100*26^(1/2) + 510
-10*10405^(1/2)
10*10405^(1/2)The Rosser matrix is not a typical example; it is rare for a full 8-by-8 matrix to have a characteristic polynomial that factors into such simple form. If you change the two "corner" elements of R from 29 to 30 with the commands
S = R; S(1,8) = 30; S(8,1) = 30;
and then try
p = poly(S)
you find
p = x^8 - 4040*x^7 + 5079941*x^6 + 82706090*x^5... - 5327831918568*x^4 + 4287832912719760*x^3... - 1082699388411166000*x^2 + 51264008540948000*x... + 40250968213600000
You also find that factor(p) is p itself. That is, the characteristic polynomial cannot be factored over the rationals.
For this modified Rosser matrix
F = eig(S)
returns
F =
1020.420188201504727818545749884
1019.9935501291629257348091808173
1019.5243552632016358324933278291
1000.1206982933841335712817075454
999.94691786044276755320289228602
0.21803980548301606860857564424981
-0.17053529728768998575200874607757
-1020.05321425589151659318942526Notice that these values are close to the eigenvalues of the original Rosser matrix. Further, the numerical values of F are a result of MuPAD software's floating-point arithmetic. Consequently, different settings of digits do not alter the number of digits to the right of the decimal place.
It is also possible to try to compute eigenvalues of symbolic matrices, but closed form solutions are rare. The Givens transformation is generated as the matrix exponential of the elementary matrix
Symbolic Math Toolbox commands
syms t A = sym([0 1; -1 0]); G = expm(t*A)
return
G =
[ (1/exp(t*i))/2 + exp(t*i)/2,
((1/exp(t*i))*i)/2 - (exp(t*i)*i)/2]
[-((1/exp(t*i))*i)/2 + (exp(t*i)*i)/2,
(1/exp(t*i))/2 + exp(t*i)/2]You can simplify this expression with the simple command:
[G,how] = simple(G)
G = [ cos(t), sin(t)] [ -sin(t), cos(t)] how = simplify
Next, the command
g = eig(G)
produces
g = cos(t) - sin(t)*i cos(t) + sin(t)*i
You can use simple to simplify this form of g:
[g,how] = simple(g)
g = 1/exp(t*i) exp(t*i) how = rewrite(exp)
The Jordan canonical form results from attempts to diagonalize a matrix by a similarity transformation. For a given matrix A, find a nonsingular matrix V, so that inv(V)*A*V, or, more succinctly, J = V\A*V, is "as close to diagonal as possible." For almost all matrices, the Jordan canonical form is the diagonal matrix of eigenvalues and the columns of the transformation matrix are the eigenvectors. This always happens if the matrix is symmetric or if it has distinct eigenvalues. Some nonsymmetric matrices with multiple eigenvalues cannot be diagonalized. The Jordan form has the eigenvalues on its diagonal, but some of the superdiagonal elements are one, instead of zero. The statement
J = jordan(A)
computes the Jordan canonical form of A. The statement
[V,J] = jordan(A)
also computes the similarity transformation. The columns of V are the generalized eigenvectors of A.
The Jordan form is extremely sensitive to perturbations. Almost any change in A causes its Jordan form to be diagonal. This makes it very difficult to compute the Jordan form reliably with floating-point arithmetic. It also implies that A must be known exactly (i.e., without roundoff error, etc.). Its elements must be integers, or ratios of small integers. In particular, the variable-precision calculation, jordan(vpa(A)), is not allowed.
For example, let
A = sym([12,32,66,116;-25,-76,-164,-294;
21,66,143,256;-6,-19,-41,-73])A = [ 12, 32, 66, 116] [ -25, -76, -164, -294] [ 21, 66, 143, 256] [ -6, -19, -41, -73]
Then
[V,J] = jordan(A)
produces
V = [ 4, -2, 4, 3] [ -6, 8, -11, -8] [ 4, -7, 10, 7] [ -1, 2, -3, -2] J = [ 1, 1, 0, 0] [ 0, 1, 0, 0] [ 0, 0, 2, 1] [ 0, 0, 0, 2]
Therefore A has a double eigenvalue at 1, with a single Jordan block, and a double eigenvalue at 2, also with a single Jordan block. The matrix has only two eigenvectors, V(:,1) and V(:,3). They satisfy
A*V(:,1) = 1*V(:,1) A*V(:,3) = 2*V(:,3)
The other two columns of V are generalized eigenvectors of grade 2. They satisfy
A*V(:,2) = 1*V(:,2) + V(:,1) A*V(:,4) = 2*V(:,4) + V(:,3)
In mathematical notation, with vj = v(:,j), the columns of V and eigenvalues satisfy the relationships
Only the variable-precision numeric computation of the complete singular vector decomposition is available in the toolbox. One reason for this is that the formulas that result from symbolic computation are usually too long and complicated to be of much use. If A is a symbolic matrix of floating-point or variable-precision numbers, then
S = svd(A)
computes the singular values of A to an accuracy determined by the current setting of digits. And
[U,S,V] = svd(A);
produces two orthogonal matrices, U and V, and a diagonal matrix, S, so that
A = U*S*V';
Consider the n-by-n matrix A with elements defined by A(i,j) = 1/(i - j + 1/2). The most obvious way of generating this matrix is
n = 5;
for i=1:n
for j=1:n
A(i,j) = sym(1/(i-j+1/2));
end
endFor n = 5, the matrix is
A
A = [ 2, -2, -2/3, -2/5, -2/7] [ 2/3, 2, -2, -2/3, -2/5] [ 2/5, 2/3, 2, -2, -2/3] [ 2/7, 2/5, 2/3, 2, -2] [ 2/9, 2/7, 2/5, 2/3, 2]
It turns out many of the singular values of these matrices are close to π.
The most efficient way to generate the matrix is
n = 5; [J,I] = meshgrid(1:n); A = sym(1./(I - J+1/2));
Since the elements of A are the ratios of small integers, vpa(A) produces a variable-precision representation, which is accurate to digits precision. Hence
S = svd(vpa(A))
computes the desired singular values to full accuracy. With n = 16 and digits(30), the result is
S = 3.14159265358979323846255035973 3.14159265358979323843066846713 3.14159265358979323325290142782 3.14159265358979270342635559052 3.1415926535897543920684990722 3.14159265358767361712392612382 3.14159265349961053143856838564 3.14159265052654880815569479613 3.14159256925492306470284863101 3.14159075458605848728982577118 3.1415575435991808369105065826 3.14106044663470063805218371923 3.13504054399744654843898901261 3.07790297231119748658424727353 2.69162158686066606774782763593 1.20968137605668985332455685355
Compare S with pi, the floating-point representation of π. In the vector below, the first element is computed by subtraction with variable-precision arithmetic and then converted to a double:
format long; double(pi*ones(16,1)-S)
The results are
ans = 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000001 0.000000000000039 0.000000000002120 0.000000000090183 0.000000003063244 0.000000084334870 0.000001899003735 0.000035109990612 0.000532206955093 0.006552109592347 0.063689681278596 0.449971066729127 1.931911277533103
Since the relative accuracy of pi is pi*eps, which is 6.9757e-16, the result confirms the suspicion that four of the singular values of the 16-by-16 example equal π to floating-point accuracy.
This example applies several numeric, symbolic, and graphic techniques to study the behavior of matrix eigenvalues as a parameter in the matrix is varied. This particular setting involves numerical analysis and perturbation theory, but the techniques illustrated are more widely applicable.
In this example, you consider a 3-by-3 matrix A whose
eigenvalues are 1, 2, 3. First, you perturb A by
another matrix E and parameter
. As t increases
from 0 to 10-6, the eigenvalues
,
,
change to
,
,
.
This, in turn, means that for some value of
, the perturbed matrix A(t) = A + tE has a double
eigenvalue
. The example shows how to find
the value of t, called τ,
where this happens.
The starting point is a MATLAB test example, known as gallery(3).
A = gallery(3)
A = -149 -50 -154 537 180 546 -27 -9 -25
This is an example of a matrix whose eigenvalues are sensitive to the effects of roundoff errors introduced during their computation. The actual computed eigenvalues may vary from one machine to another, but on a typical workstation, the statements
format long e = eig(A)
produce
e = 1.000000000010722 1.999999999991790 2.999999999997399
Of course, the example was created so that its eigenvalues are actually 1, 2, and 3. Note that three or four digits have been lost to roundoff. This can be easily verified with the toolbox. The statements
B = sym(A); e = eig(B)' p = poly(B) f = factor(p)
produce
e = [1, 2, 3] p = x^3 - 6*x^2 + 11*x - 6 f = (x - 3)*(x - 1)*(x - 2)
Are the eigenvalues sensitive to the perturbations caused by roundoff error because they are "close together"? Ordinarily, the values 1, 2, and 3 would be regarded as "well separated." But, in this case, the separation should be viewed on the scale of the original matrix. If A were replaced by A/1000, the eigenvalues, which would be .001, .002, .003, would "seem" to be closer together.
But eigenvalue sensitivity is more subtle than just "closeness." With a carefully chosen perturbation of the matrix, it is possible to make two of its eigenvalues coalesce into an actual double root that is extremely sensitive to roundoff and other errors.
One good perturbation direction can be obtained from the outer product of the left and right eigenvectors associated with the most sensitive eigenvalue. The following statement creates the perturbation matrix:
E = [130,-390,0;43,-129,0;133,-399,0]
E =
130 -390 0
43 -129 0
133 -399 0The perturbation can now be expressed in terms of a single, scalar parameter t. The statements
syms x t A = A + t*E
replace A with the symbolic representation of its perturbation:
A = [130*t - 149, - 390*t - 50, -154] [ 43*t + 537, 180 - 129*t, 546] [ 133*t - 27, - 399*t - 9, -25]
Computing the characteristic polynomial of this new A
p = simple(poly(A))
gives
p = 11*x - 1221271*t - x^2*(t + 6) + 492512*t*x + x^3 - 6
p is a cubic in x whose coefficients vary linearly with t.
It turns out that when t is varied over a very small interval, from 0 to 1.0e–6, the desired double root appears. This can best be seen graphically. The first figure shows plots of p, considered as a function of x, for three different values of t: t = 0, t = 0.5e–6, and t = 1.0e–6. For each value, the eigenvalues are computed numerically and also plotted:
x = .8:.01:3.2; for k = 0:2 c = sym2poly(subs(p,t,k*0.5e-6)); y = polyval(c,x); lambda = eig(double(subs(A,t,k*0.5e-6))); subplot(3,1,3-k) plot(x,y,'-',x,0*x,':',lambda,0*lambda,'o') axis([.8 3.2 -.5 .5]) text(2.25,.35,['t = ' num2str( k*0.5e-6 )]); end
The bottom subplot shows the unperturbed polynomial, with its three roots at 1, 2, and 3. The middle subplot shows the first two roots approaching each other. In the top subplot, these two roots have become complex and only one real root remains.
The next statements compute and display the actual eigenvalues
e = eig(A); ee = subexpr(e);
sigma = (1221271*t)/2 + (t + 6)^3/27 - ((492512*t + 11)*(t + 6))/6 +... (((492512*t)/3 - (t + 6)^2/9 + 11/3)^3 + ((1221271*t)/2 +... (t + 6)^3/27 - ((492512*t + 11)*(t + 6))/6 + 3)^2)^(1/2) + 3
pretty(ee)
showing that e(2) and e(3) form a complex conjugate pair:
+- -+ | 2 | | 1 492512 t (t + 6) 11 | | - -------- - -------- + -- | | t 3 3 9 3 | | - + sigma - ------------------------ + 2 | | 3 1 | | - | | 3 | | sigma | | | | / 2 \ | | | 1 492512 t (t + 6) 11 | | | | - -------- - -------- + -- | | | 1/2 | 3 3 9 3 | | | 3 | sigma + ------------------------ | i | | 1 2 | 1 | | | - 492512 t (t + 6) 11 | - | | | 3 -------- - -------- + -- | 3 | | | t sigma 3 9 3 \ sigma / | | - - ------ + ------------------------ + 2 - -------------------------------------------- | | 3 2 1 2 | | - | | 3 | | 2 sigma | | | | / 2 \ | | | 1 492512 t (t + 6) 11 | | | | - -------- - -------- + -- | | | 1/2 | 3 3 9 3 | | | 3 | sigma + ------------------------ | i | | 1 2 | 1 | | | - 492512 t (t + 6) 11 | - | | | 3 -------- - -------- + -- | 3 | | | t sigma 3 9 3 \ sigma / | | - - ------ + ------------------------ + 2 + -------------------------------------------- | | 3 2 1 2 | | - | | 3 | | 2 sigma | +- -+ |
Next, the symbolic representations of the three eigenvalues are evaluated at many values of t
tvals = (2:-.02:0)' * 1.e-6;
r = size(tvals,1);
c = size(e,1);
lambda = zeros(r,c);
for k = 1:c
lambda(:,k) = double(subs(e(k),t,tvals));
end
plot(lambda,tvals)
xlabel('\lambda'); ylabel('t');
title('Eigenvalue Transition')to produce a plot of their trajectories.
Above t = 0.8e-6, the graphs of two of the eigenvalues intersect, while below t = 0.8e–6, two real roots become a complex conjugate pair. What is the precise value of t that marks this transition? Let τ denote this value of t.
One way to find the exact value of τ involves polynomial discriminants. The discriminant of a quadratic polynomial is the familiar quantity under the square root sign in the quadratic formula. When it is negative, the two roots are complex.
There is no discrim function in the toolbox, but there is one in the MuPAD language. The statement
doc(symengine,'discrim')
gives the MuPAD help for the function.
This shows that the discrim function is in the polylib library. Use these commands
syms a b c x evalin(symengine,'polylib::discrim(a*x^2+b*x+c, x)')
to show the generic quadratic's discriminant, b2 - 4ac:
ans = b^2 - 4*a*c
The discriminant for the perturbed cubic characteristic polynomial is obtained, using
discrim = feval(symengine,'polylib::discrim',p,x)
which produces
discrim = 242563185060*t^4 - 477857003880091920*t^3 +... 1403772863224*t^2 - 5910096*t + 4
The quantity τ is one of the four roots of this quartic. You can find a numeric value for τ with the following code.
s = solve(discrim); tau = vpa(s)
tau =
1970031.04061804553618913725474883634597991201389
0.000000783792490596794010485879469854518820556090553664
0.00000107692481604921513807537160160597784208236311263 - 0.00000308544636502289065492747*i
0.00000308544636502289065492746538275636180217710757295*i + 0.00000107692481604921513807537160160597784249167873707 |
Of the four solutions, you know that
tau = tau(2)
is the transition point
tau = 0.00000078379249059679401048084
because it is closest to the previous estimate.
A more generally applicable method for finding τ is based on the fact that, at a double root, both the function and its derivative must vanish. This results in two polynomial equations to be solved for two unknowns. The statement
sol = solve(p,diff(p,'x'))
solves the pair of algebraic equations p = 0 and dp/dx = 0 and produces
sol =
t: [4x1 sym]
x: [4x1 sym]Find τ now by
format short tau = double(sol.t(2))
which reveals that the second element of sol.t is the desired value of τ:
tau = 7.8379e-007
Therefore, the second element of sol.x
sigma = double(sol.x(2))
is the double eigenvalue
sigma =
1.5476To verify that this value of τ does
indeed produce a double eigenvalue at
, substitute τ for t in
the perturbed matrix A(t) = A + tE and find the
eigenvalues of A(t). That is,
e = eig(double(subs(A, t, tau)))
e =
1.5476
1.5476
2.9048confirms that
is a double eigenvalue
of A(t) for t = 7.8379e–07.
 | Variable-Precision Arithmetic | Solving Equations |  |
See how symbolic computations can help you find analytical solutions to math and engineering problems.
Get free kit| © 1984-2010- The MathWorks, Inc. - Site Help - Patents - Trademarks - Privacy Policy - Preventing Piracy - RSS |