Re: Solving A x = Lam B x with A & B singular
- To: mathgroup at smc.vnet.net
- Subject: [mg15474] Re: Solving A x = Lam B x with A & B singular
- From: Daniel Lichtblau <danl>
- Date: Mon, 18 Jan 1999 04:22:14 -0500
- Organization: Wolfram Research, Inc.
- References: <77f0qi$m1u@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Below are a few crude ways that may or may not work for relatively simple examples. I am very doubtful that anything less that full-blown generalized eigenvalues codes, a la Lapack, will be effective for large and/or numerically dicey problems. First method: set up an appropriate interpolating polynomial and find the roots. genEigvals1[A_, B_, var_] := Module[{n=Length[A], vals, gendet, rank}, rank = Length[SingularValues[B][[2]]]; vals = Map[Det[A-#*B]&, Table[j, {j, 1, rank+1}]]; gendet = InterpolatingPolynomial[vals, var]; Chop[Apply[List,Map[Last,NRoots[gendet==0, var]]]] ] Second method: Perturb nearby nonsingular examples and then average them. genEigvals2[A_, B_, var_] := Module[ {n=Length[A], dd, eigvals1, eigvals2, eigvals3, eigvals4}, dd = 10^(-3)*IdentityMatrix[n]; eigvals1 = Eigenvalues[Inverse[B+dd].A]; eigvals2 = Eigenvalues[Inverse[B-dd].A]; eigvals3 = Eigenvalues[Inverse[B+I*dd].A]; eigvals4 = Eigenvalues[Inverse[B-I*dd].A]; Chop[(eigvals1 + eigvals2 + eigvals3 + eigvals4)/4] ] Third method: Set up and explicitly solve the generalized eigen equation. genEigvals3[A_, B_, var_] := Module[ {n=Length[A], j, k, coeffs, c, len, vec, lam, tmpvec}, coeffs = Array[c, n-1]; For [j=1, j<=n, j++, vec[j] = Table [If [k==j, 1, 0], {k,n}]]; tmpvec = vec[n] + Sum[c[j]*vec[j], {j,n-1}]; Chop[lam /. Solve[A . tmpvec == lam*B.tmpvec, Prepend[coeffs,lam]]] ] For a simple example, I'll use a pair of singular 3x3 matrices. A = N[{{1,2,3}, {4,5,6}, {7,8,9}}]; B = {{-1, 6, 7}, {8, 2, -5.}, {7, 8, 2}}; In[15]:= gvals1 = genEigvals1[A, B, x] Out[15]= {-0.891892, 0} In[16]:= gvals2 = genEigvals2[A,B,x] Out[16]= {2.68044, -0.891892, 0} In[17]:= gvals3 = genEigvals3[A,B,x] Out[17]= {-0.891892, 0} Note that the second method will give junk even for this example because moving to a nonsingular example forces a third eigenvalue. Also note that it is assuming the standard sorting will keep corresponding eigenvalues in the same position, which might fail if there is a pair of generalized eigenvalues with nearly equal real parts (if I correctly remember our ordering of eigenvalues). The third method will fail if there are eigenvectors with last component equal to zero. I am fairly certain there are other ways in which these may give wrong or incomplete results. Daniel Lichtblau Wolfram Research