Bill:
Your matrix is wrong. Unless I am misreading the paper, the first matrix is
m=:4 4 $ 1,(%3),5 1 3 1 5 1, (%5), (%5), 1, (%5), 1 1 5 1
not
n=:4 4 $ 1,(%3),5 1 3 1 5 1, (%2), (%2), 1, (%2), 1 1 5 1
which is what you have got.
Then, where power is the code I posted before and la is the same from LAPACK,
la=:3 : '{."1 (2{::dgeev_jlapack_ y.)'
I get
(%+/) power m
0.232323 0.409949 0.0596214 0.298107
(%+/) la m
0.232323 0.409949 0.0596214 0.298107
which is close but not equal to Coyle's result. I still do not quite
undertand what "standard methods" he is using to calculate eigenvalues.
The power method is the easiest for a dominant eigenvalue (if it works).
LAPACK uses a completely different method (the matrix is balanced and put
in upper Hessenberg form). These (standard) methods give identical
results which are different from Coyle's.
Please feel free to use any code I have posted for any purpose.
Best wishes,
John
Bill Harris wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> "John Randall" <[EMAIL PROTECTED]> writes:
>
>> Here is a slightly more comprehensible form of the power method code I
>> posted yesterday. It will correctly calculate the principal
>> eigenvectors
>> in Coyle's paper.
>
> Actually, I get the same answers on the 3x3 matrices but not the
> original 4x4.
>
> cpm
> 1 0.33333333 5 1
> 3 1 5 1
> 0.5 0.5 1 0.5
> 1 1 5 1
> (%+/) power cpm
> 0.23354981 0.37461869 0.1091339 0.28269759
>
> Coyle's eigenvector: (0.232, 0.402, 0.061, 0.305)
>
> It's not far off, though.
>
> Bill
> - --
> Bill Harris http://facilitatedsystems.com/weblog/
> Facilitated Systems Everett, WA 98208 USA
> http://facilitatedsystems.com/ phone: +1 425 337-5541
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.1 (MingW32)
> Comment: For more information, see http://www.gnupg.org
>
> iD8DBQFEB1eH3J3HaQTDvd8RAosdAJ4mJOmbK2+xpqOWuYoWvaJ4M5q0rQCeLUhr
> ryGA+X7CVdk/thvH0GW1h48=
> =sR+1
> -----END PGP SIGNATURE-----
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
