## File: Eigenvalue-and-Eigenvector-Examples.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (292 lines) | stat: -rw-r--r-- 9,512 bytes parent folder | download
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292` `````` GNU Scientific Library – Reference Manual: Eigenvalue and Eigenvector Examples

15.8 Examples

The following program computes the eigenvalues and eigenvectors of the 4-th order Hilbert matrix, H(i,j) = 1/(i + j + 1).

#include <stdio.h> #include <gsl/gsl_math.h> #include <gsl/gsl_eigen.h>  int main (void) {   double data[] = { 1.0  , 1/2.0, 1/3.0, 1/4.0,                     1/2.0, 1/3.0, 1/4.0, 1/5.0,                     1/3.0, 1/4.0, 1/5.0, 1/6.0,                     1/4.0, 1/5.0, 1/6.0, 1/7.0 };    gsl_matrix_view m      = gsl_matrix_view_array (data, 4, 4);    gsl_vector *eval = gsl_vector_alloc (4);   gsl_matrix *evec = gsl_matrix_alloc (4, 4);    gsl_eigen_symmv_workspace * w =      gsl_eigen_symmv_alloc (4);      gsl_eigen_symmv (&m.matrix, eval, evec, w);    gsl_eigen_symmv_free (w);    gsl_eigen_symmv_sort (eval, evec,                          GSL_EIGEN_SORT_ABS_ASC);      {     int i;      for (i = 0; i < 4; i++)       {         double eval_i             = gsl_vector_get (eval, i);         gsl_vector_view evec_i             = gsl_matrix_column (evec, i);          printf ("eigenvalue = %g\n", eval_i);         printf ("eigenvector = \n");         gsl_vector_fprintf (stdout,                              &evec_i.vector, "%g");       }   }    gsl_vector_free (eval);   gsl_matrix_free (evec);    return 0; }

Here is the beginning of the output from the program,

\$ ./a.out  eigenvalue = 9.67023e-05 eigenvector =  -0.0291933 0.328712 -0.791411 0.514553 ...

This can be compared with the corresponding output from GNU OCTAVE,

octave> [v,d] = eig(hilb(4)); octave> diag(d)   ans =     9.6702e-05    6.7383e-03    1.6914e-01    1.5002e+00  octave> v  v =     0.029193   0.179186  -0.582076   0.792608   -0.328712  -0.741918   0.370502   0.451923    0.791411   0.100228   0.509579   0.322416   -0.514553   0.638283   0.514048   0.252161

Note that the eigenvectors can differ by a change of sign, since the sign of an eigenvector is arbitrary.

The following program illustrates the use of the nonsymmetric eigensolver, by computing the eigenvalues and eigenvectors of the Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4).

#include <stdio.h> #include <gsl/gsl_math.h> #include <gsl/gsl_eigen.h>  int main (void) {   double data[] = { -1.0, 1.0, -1.0, 1.0,                     -8.0, 4.0, -2.0, 1.0,                     27.0, 9.0, 3.0, 1.0,                     64.0, 16.0, 4.0, 1.0 };    gsl_matrix_view m      = gsl_matrix_view_array (data, 4, 4);    gsl_vector_complex *eval = gsl_vector_complex_alloc (4);   gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);    gsl_eigen_nonsymmv_workspace * w =      gsl_eigen_nonsymmv_alloc (4);      gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);    gsl_eigen_nonsymmv_free (w);    gsl_eigen_nonsymmv_sort (eval, evec,                             GSL_EIGEN_SORT_ABS_DESC);      {     int i, j;      for (i = 0; i < 4; i++)       {         gsl_complex eval_i             = gsl_vector_complex_get (eval, i);         gsl_vector_complex_view evec_i             = gsl_matrix_complex_column (evec, i);          printf ("eigenvalue = %g + %gi\n",                 GSL_REAL(eval_i), GSL_IMAG(eval_i));         printf ("eigenvector = \n");         for (j = 0; j < 4; ++j)           {             gsl_complex z =                gsl_vector_complex_get(&evec_i.vector, j);             printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));           }       }   }    gsl_vector_complex_free(eval);   gsl_matrix_complex_free(evec);    return 0; }

Here is the beginning of the output from the program,

\$ ./a.out  eigenvalue = -6.41391 + 0i eigenvector =  -0.0998822 + 0i -0.111251 + 0i 0.292501 + 0i 0.944505 + 0i eigenvalue = 5.54555 + 3.08545i eigenvector =  -0.043487 + -0.0076308i 0.0642377 + -0.142127i -0.515253 + 0.0405118i -0.840592 + -0.00148565i ...

This can be compared with the corresponding output from GNU OCTAVE,

octave> [v,d] = eig(vander([-1 -2 3 4])); octave> diag(d) ans =    -6.4139 + 0.0000i    5.5456 + 3.0854i    5.5456 - 3.0854i    2.3228 + 0.0000i  octave> v v =   Columns 1 through 3:    -0.09988 + 0.00000i  -0.04350 - 0.00755i  -0.04350 + 0.00755i   -0.11125 + 0.00000i   0.06399 - 0.14224i   0.06399 + 0.14224i    0.29250 + 0.00000i  -0.51518 + 0.04142i  -0.51518 - 0.04142i    0.94451 + 0.00000i  -0.84059 + 0.00000i  -0.84059 - 0.00000i   Column 4:    -0.14493 + 0.00000i    0.35660 + 0.00000i    0.91937 + 0.00000i    0.08118 + 0.00000i

Note that the eigenvectors corresponding to the eigenvalue 5.54555 + 3.08545i differ by the multiplicative constant 0.9999984 + 0.0017674i which is an arbitrary phase factor of magnitude 1.

``````