Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Polynomial eigenvalue problem

• To: mathgroup at smc.vnet.net
• Subject: [mg58725] Re: [mg58699] Polynomial eigenvalue problem
• From: "Carl K. Woll" <carl at woll2woll.com>
• Date: Fri, 15 Jul 2005 14:14:07 -0400 (EDT)
• References: <200507150702.DAA23070@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

----- Original Message ----- 
From: "Alan" <info at optioncity.REMOVETHIS.net>
To: mathgroup at smc.vnet.net
Subject: [mg58725] [mg58699] Polynomial eigenvalue problem


> Does anyone know a Mathematica solution
> for the polynomial eigenvalue problem (PEP)?
>
> The PEP problem of order p . N is:
> Find the eigenvalues lambda, and eigenvectors v (each an N-vector) that
> solve
>
> (A_0 + lambda A_1 + lambda^2 A_2 + ... + lambda^p A_p) x = 0,
>
> where:
> p is a given integer,
> each coefficients A_i is a given N x N (real or complex) square matrix.
>
> There are p . N solution pairs (lambda, x) if the problem is not 
> ill-posed.
> If both A_0 and A_p are singular the problem is potentially ill-posed.
>
> Thanks for any help,
> alan
>

In version 5, the linear problem can be directly solved using Eigensystem if 
A_0 is nonsingular:

Eigensystem[{A_0, -A_1}]

For higher order polynomials, one idea is to proceed the same way as one 
does with the usual eigenvalue problem. First some test data:

randmatrix[n_] := Table[Random[], {n}, {n}]

SeedRandom[1]
a[0] = randmatrix[3];
a[1] = randmatrix[3];
a[2] = randmatrix[3];
mat = a[0] + lambda a[1] + lambda^2 a[2];

Get the eigenvalues by setting determinant to zero:

eigs = Sort[NSolve[Det[mat]]];

Get corresponding eigenvectors by using NullSpace:

eigvecs = Join @@ NullSpace /@ (mat /. Union[eigs]);

The Sort and Union commands are there in case there are repeated 
eigenvalues.

Let's check that our eigenvalues and eigenvectors satisfy your equation:

In[21]:=
Table[(mat/.eigs[[i]]) . eigvecs[[i]], {i,Length[eigs]}]//Chop

Out[21]=
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

Let's compare our results for the linear case:

mat = a[0] + lambda a[1];

In[23]:=
NSolve[Det[mat]]

Out[23]=
{{lambda -> -1.29964}, {lambda -> -0.97271}, {lambda -> 1.38013}}

In[24]:=
Eigenvalues[{a[0],-a[1]}]

Out[24]=
{1.38013, -1.29964, -0.97271}

Carl Woll
Wolfram Research 

• References:
  • Polynomial eigenvalue problem
    • From: "Alan" <info@optioncity.REMOVETHIS.net>