How to calculate the Eigenvalue of a symbolic matrix

Question

I have a problem to calculate eigenvalues of a symbolic matrix. here is the code is matlab format:

clc
clear all

i=1.824e-6;
j=0.11e-6;
c=7.342e-9;
e=200e9;
g=48.8e9;
a=0.075;
l=1;

syms p

k=[12*e*i/l^3,-6*e*i/l^2,0,0,-12*e*i/l^3,-6*e*i/l^2,0,0;-6*e*i/l^2,4*e*i/l,0,0,6*e*i/l^2,2*e*i/l,0,0;0,0,1.2*g*j/l+12*e*c/l^3,-0.1*g*j-6*e*c/l^2,0,0,-1.2*g*j/l-12*e*c/l^3,-1.2*g*j/l-12*e*c/l^3;0,0,-0.1*g*j-6*e*c/l^2,2/15*g*j*l+4*e*c/l,0,0,0.1*g*j+6*e*c/l^2,-1/30*g*j*l+2*e*c/l;-12*e*i/l^3,6*e*i/l^2,0,0,12*e*i/l^3,6*e*i/l^2,0,0;-6*e*i/l^2,2*e*i/l,0,0,6*e*i/l^2,4*e*i/l,0,0;0,0,-1.2*g*j/l-12*e*c/l^3,0.1*g*j+6*e*c/l^2,0,0,1.2*g*j/l+12*e*c/l^3,0.1*g*j+6*e*c/l^2;0,0,-0.1*g*j-6*e*c/l^2,-1/30*g*j*l+2*e*c/l,0,0,0.1*g*j+6*e*c/l^2,2/15*g*j*l+4*e*c/l];

kG=p/840*[0,0,198*l,-15*l^2,0,0,-198*l,-15*l^2;0,0,-169*l^2,18*l^3,0,0,29*l^2,-3*l^3;198*l,-169*l^2,-312*a*l,44*a*l^2,-198*l,-29*l^2,-108*a*l,-26*a*l^2;-15*l^2,18*l^3,44*a*l^2,-8*a*l^3,15*l^2,-3*l^3,26*a*l^2,6*a*l^3;0,0,-198*l,15*l^2,0,0,198*l,15*l^2;0,0,-29*l^2,-3*l^3,0,0,169*l^2,18*l^3;-198*l,29*l^2,-108*a*l,26*a*l^2,198*l,169*l^2,-312*a*l,-44*a*l^2;-15*l^2,-3*l^3,-26*a*l^2,6*a*l^3,15*l^2,18*l^3,-44*a*l^2,-8*a*l^3];

eig(k+kG)

I want a numerical result. matlab gives me as "root(...)". I heard maple can solve this problem but I am amateur on maple. please help me.

Answer 1

If you convert all your floating-point coefficients to exact rationals then it appears that k+kG has eigenvalue zero with multiplicity of two.

But the characteristic polynomial that remains (after dividing out lambda^2) is of degree 6, and that does not appear to factor directly.

So you will need to substitute numeric values for p, in order to obtain the remaining eigenvalues as floating-point. Perhaps the following code will help get you started.

restart;
Digits := 15:
par := [i=1.824e-6, j=0.11e-6, c=7.342e-9,
        e=200e9, g=48.8e9, a=0.075, l=1]:
epar := convert(par, rational, exact);

k := Matrix([[12*e*i/l^3,-6*e*i/l^2,0,0,-12*e*i/l^3,-6*e*i/l^2,0,0],
            [-6*e*i/l^2,4*e*i/l,0,0,6*e*i/l^2,2*e*i/l,0,0],
            [0,0,1.2*g*j/l+12*e*c/l^3,-0.1*g*j-6*e*c/l^2,0,0,-1.2*g*j/l-12*e*c/l^3,-1.2*g*j/l-12*e*c/l^3],
            [0,0,-0.1*g*j-6*e*c/l^2,2/15*g*j*l+4*e*c/l,0,0,0.1*g*j+6*e*c/l^2,-1/30*g*j*l+2*e*c/l],
            [-12*e*i/l^3,6*e*i/l^2,0,0,12*e*i/l^3,6*e*i/l^2,0,0],
            [-6*e*i/l^2,2*e*i/l,0,0,6*e*i/l^2,4*e*i/l,0,0],
            [0,0,-1.2*g*j/l-12*e*c/l^3,0.1*g*j+6*e*c/l^2,0,0,1.2*g*j/l+12*e*c/l^3,0.1*g*j+6*e*c/l^2],
            [0,0,-0.1*g*j-6*e*c/l^2,-1/30*g*j*l+2*e*c/l,0,0,0.1*g*j+6*e*c/l^2,2/15*g*j*l+4*e*c/l]]):

k := map(convert,k,rational,exact):

kG := p/840*Matrix([[0,0,198*l,-15*l^2,0,0,-198*l,-15*l^2],
                    [0,0,-169*l^2,18*l^3,0,0,29*l^2,-3*l^3],
                    [198*l,-169*l^2,-312*a*l,44*a*l^2,-198*l,-29*l^2,-108*a*l,-26*a*l^2],
                    [-15*l^2,18*l^3,44*a*l^2,-8*a*l^3,15*l^2,-3*l^3,26*a*l^2,6*a*l^3],
                    [0,0,-198*l,15*l^2,0,0,198*l,15*l^2],
                    [0,0,-29*l^2,-3*l^3,0,0,169*l^2,18*l^3],
                    [-198*l,29*l^2,-108*a*l,26*a*l^2,198*l,169*l^2,-312*a*l,-44*a*l^2],
                    [-15*l^2,-3*l^3,-26*a*l^2,6*a*l^3,15*l^2,18*l^3,-44*a*l^2,-8*a*l^3]]):

M := k+kG:
Q := eval(M, epar):

with(LinearAlgebra):

cp := CharacteristicPolynomial(Q, lambda):
cp := lambda^2*collect(cp/lambda^2,lambda,factor);

solve(cp);

sort([fsolve(eval(cp, p=2.0), complex)], (a,b)->abs(a)>abs(b));
sort(Eigenvalues(eval(Q, p=2.0)), (a,b)->abs(a)>abs(b));

sort([fsolve(eval(cp, p=3000.0), complex)], (a,b)->Re(a)>Re(b));
sort(Eigenvalues(eval(Q, p=3000.0)), (a,b)->Re(a)>Re(b));

I could mention that kG is symmetric, and k is almost symmetric. The entries k[3,8] and k[8,3] differ. Is that intentional, or a mistake? (If k is intended to be symmetric then please clarify, because it means that k+kG can be cast as a Maple Matrix with shape=symmetric and the floating-point eigenvalue computation will be a little faster, produce no 0.0*I imaginary artefacts, and return sorted by default.)

[edit] Following clarification by the Original Poster I make an adjustment below to replace the former definition of the [3,8] entry of Matrix k by the value for its [8,3] entry.

That makes Matrix k symmetric.

That also allows the characteristic polynomial of the exact rational form of k+kG to factor symbolically w.r.t. lambda. This results in 8 explicit formulas (in terms of parameter p) for the eigenvalues. Note however that these involve radicals and the imaginary unit, and immediate substitution of a floating-point value for parameter p results in small nonzero imaginary component in the results (as numeric roundoff artefact).

It thus still may be superior to use a dedicated root-finder on the characteristic polynomial, or compute the eigenvalues numerically, following substitution at floating-point values of parameter p. These need not produce nonzero imaginary component artefacts as noise due to roundoff error.

Here is revised code, with the now symmetric Matrix k.

restart;
Digits := 15:
par := [i=1.824e-6, j=0.11e-6, c=7.342e-9,
        e=200e9, g=48.8e9, a=0.075, l=1]:
epar := convert(par, rational, exact):

k := Matrix([[12*e*i/l^3,-6*e*i/l^2,0,0,-12*e*i/l^3,-6*e*i/l^2,0,0],
            [-6*e*i/l^2,4*e*i/l,0,0,6*e*i/l^2,2*e*i/l,0,0],
            [0,0,1.2*g*j/l+12*e*c/l^3,-0.1*g*j-6*e*c/l^2,0,0,-1.2*g*j/l-12*e*c/l^3,-0.1*g*j-6*e*c/l^2],
            [0,0,-0.1*g*j-6*e*c/l^2,2/15*g*j*l+4*e*c/l,0,0,0.1*g*j+6*e*c/l^2,-1/30*g*j*l+2*e*c/l],
            [-12*e*i/l^3,6*e*i/l^2,0,0,12*e*i/l^3,6*e*i/l^2,0,0],
            [-6*e*i/l^2,2*e*i/l,0,0,6*e*i/l^2,4*e*i/l,0,0],
            [0,0,-1.2*g*j/l-12*e*c/l^3,0.1*g*j+6*e*c/l^2,0,0,1.2*g*j/l+12*e*c/l^3,0.1*g*j+6*e*c/l^2],
            [0,0,-0.1*g*j-6*e*c/l^2,-1/30*g*j*l+2*e*c/l,0,0,0.1*g*j+6*e*c/l^2,2/15*g*j*l+4*e*c/l]]):

k := map(convert,k,rational,exact):

kG := p/840*Matrix([[0,0,198*l,-15*l^2,0,0,-198*l,-15*l^2],
                    [0,0,-169*l^2,18*l^3,0,0,29*l^2,-3*l^3],
                    [198*l,-169*l^2,-312*a*l,44*a*l^2,-198*l,-29*l^2,-108*a*l,-26*a*l^2],
                    [-15*l^2,18*l^3,44*a*l^2,-8*a*l^3,15*l^2,-3*l^3,26*a*l^2,6*a*l^3],
                    [0,0,-198*l,15*l^2,0,0,198*l,15*l^2],
                    [0,0,-29*l^2,-3*l^3,0,0,169*l^2,18*l^3],
                    [-198*l,29*l^2,-108*a*l,26*a*l^2,198*l,169*l^2,-312*a*l,-44*a*l^2],
                    [-15*l^2,-3*l^3,-26*a*l^2,6*a*l^3,15*l^2,18*l^3,-44*a*l^2,-8*a*l^3]]):

M := k+kG:
Q := eval(M, epar):

with(LinearAlgebra):

cp := CharacteristicPolynomial(Q, lambda):
cp := collect(cp/lambda^2,lambda,factor):

# This produces R, as explicit symbolic formulas
# for the eigenvalues, in terms of parameter p
R := map[2](eval,lambda,[solve(cp)]):

# Extra precision is taken, to try and deal with
# small imaginary artefacts produced when evaluating
# those formula at specific numeric value of p.
sort(evalf[100](eval(R, p=2.0)), (a,b)->abs(a)>abs(b)):
simplify(fnormal(evalf(%)),zero);

sort([fsolve(eval(cp, p=2.0), complex)], (a,b)->abs(a)>abs(b));

sort(Eigenvalues(Matrix(eval(Q, p=2.0),shape=symmetric)), (a,b)->abs(a)>abs(b));