Author 
Message 
rjfmv science forum beginner
Joined: 30 Jan 2006
Posts: 13

Posted: Wed Jul 12, 2006 4:50 pm Post subject:
Quadratic eigenvalue problem



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(BBt);
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; 

Back to top 


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

Posted: Wed Jul 12, 2006 5:47 pm Post subject:
Re: Quadratic eigenvalue problem



In article <1152723021.483424.257370@m73g2000cwd.googlegroups.com>,
rjfmv <rjfm.vieira@gmail.com> wrote:
Quote:  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

....
Quote:  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 suspect the basic problem is that the nonzero eigenvalues are
somewhat complicated. It may be better to do something like this:
Quote:  P:= select(has, factor(lhs(eq_cct)), lambda);

12 2 2 2
P := lambda (2 lambda E~ k~ b~ + 4 b~ k~ lambda E~ h~
4 4 2 2
+ lambda G~ E~ b~ + lambda G~ E~ h~ + 2 k~ h~ lambda E~
 96 k~ G~)
Quote:  R:= RootOf(op(2,P), lambda);

4
R := RootOf((G~ E~ b~ + G~ E~ h~) _Z
2 2 2
+ (2 E~ k~ b~ + 2 k~ h~ E~ + 4 b~ k~ E~ h~) _Z  96 k~ G~)
Quote:  subs(R=lambda,nullspace(eval(Q,lambda=R)));

[
{[1, 1, 1, 1, lambda
2 2 2 2
(2 k~ b~  lambda G~ b~  lambda G~ h~  2 k~ h~ ) E~/(48
k~ G~ b~), E~
2 2 2 2
(2 k~ b~  2 k~ h~ + lambda G~ b~ + lambda G~ h~) lambda/
(48 k~ G~ h~),  lambda
2 2 2 2
(2 k~ b~  lambda G~ b~  lambda G~ h~  2 k~ h~ ) E~/(48
k~ G~ b~),  E~
2 2 2 2
(2 k~ b~  2 k~ h~ + lambda G~ b~ + lambda G~ h~) lambda/
]
(48 k~ G~ h~)]}
Robert Israel israel@math.ubc.ca
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia Vancouver, BC, Canada 

Back to top 


Google


Back to top 



The time now is Wed Nov 22, 2017 10:10 pm  All times are GMT

