Quadratic eigenvalue problem

rjfmv · science forum beginner Joined: 30 Jan 2006 Posts: 13

Dear group,

I am using Maple 7 to solve a quadratic eigenvalue problem of the type:

A*lambda^2+B*lambda+C=0

A and C are symmetric matrices; B is a skew symmetric matrix.

The matrices A,B and C are written in terms of constants: E, G, b, h, k
which are positive. I have used the command assume in order to
implement this assumption:

assume(k>0);assume(E>0);assume(G>0);assume(b>0);assume(h>0);

I have obtained the eigenvalues through characteristic equation as
follows:

det(A*lambda^2+B*lambda+C)=0

In what concerns eigenvectors two distinct situations arise:

1) zero eigenvalues

The eigenvectors associated with the zero eigenvalue are easy to obtain
through the computation of the nullspace of C. However the zero
eigenvalue has an algebraic multiplicity of 12 and the kernel of C is
made of 4 eigenvectors, therefore a Jordan chain must be computed as
follows

C*d0=0
B*d0+C*d1=0
A*d2+B*d1+C*d0=0
A*d3+B*d2+C*d1=0

for d0 <> 0 and where d0,d1,d2 and d3 are the corresponding
eigenvectors.

Since it requires some iteration, I have found difficulties in
implementing the computation of the Jordan chains symbolically. If
anyone has an idea how to solve this issue, I would be grateful for an
answer.

2) non zero eigenvalues

For non zero eigenvalues, I have tried unsuccessfully to obtain the
corresponding eigenvectors
through the solution of a linear system of equations with linsolve as
follows,

(A*lambda^2+B*lambda+C)V=0, where the only unkown is v the eigenvector.

Instead of linsolve, I have also used the command nullspace to compute
V.

However, my maple worksheet didn't solve the linear system nor the
nullspace. There wasn't any error/warning message, it just keep working
for hours without reaching a solution

Is it possible that maple is considering the symbolic constants E, G,
b, h, k as variables? and if so how can I tell maple to consider them
as constants?

I have enclosed my maple worsheet for better understanding.

Hoping that you could help me out, I thank you for your attention

Ricardo

---------------------------------------------------------------------------------------------------------------------------------

assume(k>0);assume(E>0);assume(G>0);assume(b>0);assume(h>0);assume(nu>0);

a11:=-E*(b+h)/3:
a12:=-E*b/6:
a13:=0:
a14:=-E*h/6:

a22:=-E*(b+h)/3:
a23:=-E*h/6:
a24:=0:

a33:=-E*(b+h)/3:
a34:=-E*b/6:

a44:=-E*(b+h)/3:

a55:=-G*b:
a66:=-G*h:
a77:=a55:
a88:=a66:

A:=matrix(8,8,([[a11,a12,a13,a14,0,0,0,0],
[a12,a22,a23,a24,0,0,0,0],
[a13,a23,a33,a34,0,0,0,0],
[a14,a24,a34,a44,0,0,0,0],
[0,0,0,0,a55,0,0,0],
[0,0,0,0,0,a66,0,0],
[0,0,0,0,0,0,a77,0],
[0,0,0,0,0,0,0,a88]]));

b11:=G:

b12:=-G:

b22:=G:

b23:=-G:

b33:=-G:

b34:=G:

b41:=G:

b44:=-G:

B:=matrix(8,8,([[0,0,0,0,-b11,-b12,0,0],
[0,0,0,0,0,-b22,-b23,0],
[0,0,0,0,0,0,-b33,-b34],
[0,0,0,0,-b41,0,0,-b44],
[b11,b12,0,0,0,0,0,0],
[0,b22,b23,0,0,0,0,0],
[0,0,b33,b34,0,0,0,0],
[b41,0,0,b44,0,0,0,0]])):

B:=matrix(8,8,([[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[b11,b12,0,0,0,0,0,0],
[0,b22,b23,0,0,0,0,0],
[0,0,b33,b34,0,0,0,0],
[b41,0,0,b44,0,0,0,0]])):

Bt:=transpose(B):
B:=evalm(B-Bt);

c11:=G*(1/b+1/h):

c12:=-G*(1/b):

c13:=0:

c14:=-G*(1/h):

c22:=G*(1/b+1/h):

c23:=-G*(1/h):

c24:=0:

c33:=G*(1/b+1/h):

c34:=-G*(1/b):

c44:=G*(1/b+1/h):

C:=matrix(8,8,([[c11,c12,c13,c14,0,0,0,0],
[c12,c22,c23,c24,0,0,0,0],
[c13,c23,c33,c34,0,0,0,0],
[c14,c24,c34,c44,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0]])):

C1:=matrix(8,8,([[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0],
[0,0,0,0,-b^2,b*h,b^2,-b*h],
[0,0,0,0,b*h,-h^2,-b*h,h^2],
[0,0,0,0,b^2,-b*h,-b^2,b*h],
[0,0,0,0,-b*h,h^2,b*h,-h^2]])):

C1:=evalm( k*C1):

C:=evalm(C+C1);

Q:=evalm(lambda^2*A+lambda*B+C);

eq_cct:=det(Q)=0;
sol:=solvefor[lambda](eq_cct):

N:=nops(sol);
val_pp:=matrix(N,1):

for i from 1 to N do

val_pp[i,1]:=simplify(rhs(op(i,sol)));

end do;

for i from 1 to N do

Q11[i]:=subs(lambda=val_pp[i,1],op(Q)):
z:=vector([0,0,0,0,0,0,0,0]):

d[i]:=convert(linsolve(Q11[i],z),matrix):
end do;

for i from 1 to N do

evalm(d[i]);

map(limit,evalm(d[i]),k=infinity);

end do;

Robert B. Israel · science forum Guru Joined: 24 Mar 2005 Posts: 2151

In article <1152723021.483424.257370@m73g2000cwd.googlegroups.com>,
rjfmv <rjfm.vieira@gmail.com> wrote:

Google