Iterative Estimates for Eigenvalues

Power Method
Lay Example 5.8.2

We illustrate the use of the Power Method to approximate a specific eigenvalue and eigenvector of a given matrix.
Initial data.

In[1]:=

Clear[a, x, μ] ; 

a = ( {{6, 5}, {1, 2}} ) ; 

x[0] = ({{0}, {1}}) ;

Use the Power Method to compute the μs and xs.

In[4]:=

termWithLargestAbsoluteValue[xs_] := If[Abs[Max[xs]] >Abs[Min[xs]], Max[xs], Min[xs]] ; 

n = 5 ; 

Do[μ[k] = termWithLargestAbsoluteValue[a . x[k]] ;      x[k + 1] = (1/μ[k]) a . x[k],     {k, 0, n}]

Assemble Lay's Table 2, p366.

In[7]:=

table = Transpose[Table[{x[k]//MatrixForm, a . x[k]//MatrixForm, μ[k]},  & ... ;              {k, 0, n}]//N] ;

ks = Prepend[Range[1, n], "k=0"] ; 

powerMethodTable = TableForm[table,      TableHeadings {{"x[k]", "a.x[k]", "μ[k]"}, ks}]

Out[9]//TableForm=

k=0 1 2 3 4 5
x[k] ( {{0.}, {1.}} ) ( {{1.}, {0.4}} ) ( {{1.}, {0.225}} ) ( {{1.}, {0.203509}} ) ( {{1.}, {0.2005}} ) ( {{1.}, {0.200071}} )
a.x[k] ( {{5.}, {2.}} ) ( {{8.}, {1.8}} ) ( {{7.125}, {1.45}} ) ( {{7.01754}, {1.40702}} ) ( {{7.0025}, {1.401}} ) ( {{7.00036}, {1.40014}} )
μ[k] 5.` 8.` 7.125` 7.017543859649122` 7.0025` 7.0003570153516606`

Check the results obtained from the Power Method.

In[10]:=

v1 = x[n]//Flatten//N ;

λ1 = termWithLargestAbsoluteValue[a . v1]//N ; 

a . v1 - λ1 v1

Out[12]=

{0., -0.000428444}

Check using Mathematica's Eigensystem command.

In[13]:=

{evals, evecs} = Eigensystem[a]

Out[13]=

{{7, 1}, {{5, 1}, {-1, 1}}}

Inverse Power Method
Lay Example 5.8.3

We illustrate the use of the Inverse Power Method to approximate a specific eigenvalue and eigenvector of a given matrix.
Initial data.

In[14]:=

Clear[a, x, μ, ν] ; 

a = ({{10, -8, -4}, {-8, 13, 4}, {-4, 5, 4}}) ; 

x[0] = ({{1}, {1}, {1}}) ; 

α = 1.9 ;

Use the Inverse Power Method to compute the μs, νs, and xs.
This naive implementation of the Inverse Power Method uses the inverse of the matrix (a-α id) to compute y[k].
A more efficient version would solve a matrix equation for y[k].

In[18]:=

id = IdentityMatrix[3] ;

inv = Inverse[a - α id] ; 

termWithLargestAbsoluteValue[xs_] := If[Abs[Max[xs]] >Abs[Min[xs]], Max[xs], Min[xs]] ; 

n = 6 ; 

Do[y[k] = inv . x[k] ;      μ[k] = termWithLargestAbso ... bsp;    x[k + 1] = (1/μ[k]) y[k],     {k, 0, n}]

Assemble Lay's Table 3, p368.

In[23]:=

table = Transpose[Table[{x[k]//MatrixForm, y[k]//MatrixForm, μ[k], ν[k]},  ... ;              {k, 0, n}]//N] ;

ks = Prepend[Range[1, n], "k=0"] ; 

inversePowerMethodTable = TableForm[table, TableHeadings {{"x[k]", "y[k]", "μ[k]", "ν[k]"}, ks}]

Out[25]//TableForm=

k=0 1 2 3 4 5 6
x[k] ( {{1.}, {1.}, {1.}} ) ( {{0.57359}, {0.0646492}, {1.}} ) ( {{0.505364}, {0.00445297}, {1.}} ) ( {{0.500378}, {0.000312808}, {1.}} ) ( {{0.500027}, {0.0000220072}, {1.}} ) ( {{0.500002}, {1.54846*10^-6}, {1.}} ) ( {{0.5}, {1.08953*10^-7}, {1.}} )
y[k] ( {{4.45037}, {0.501601}, {7.7588}} ) ( {{5.01306}, {0.0441721}, {9.9197}} ) ( {{5.00125}, {0.00312649}, {9.99493}} ) ( {{5.00009}, {0.000220064}, {9.99965}} ) ( {{5.00001}, {0.0000154845}, {9.99998}} ) ( {{5.}, {1.08953*10^-6}, {10.}} ) ( {{5.}, {7.66613*10^-8}, {10.}} )
μ[k] 7.758804695837776` 9.91970041059954` 9.994930860390895` 9.999646313779577` 9.999975129049835` 9.999998250104962` 9.999999876874135`
ν[k] 2.02888583218707` 2.0008094961145666` 2.0000507171053` 2.0000035369873026` 2.00000024871012` 2.0000000174989534` 2.0000000012312587`

Check the results obtained from the Inverse Power Method.

In[26]:=

v1 = x[n]//Flatten//N ;

λ1 = termWithLargestAbsoluteValue[a . v1]//N ; 

a . v1 - λ1 v1

Out[28]=

{1.74158*10^-7, 1.43951*10^-7, 0.}

Check using Mathematica's Eigensystem command.

In[29]:=

{evals, evecs} = Eigensystem[a]//N

Out[29]=

{{21.6788, 3.32122, 2.}, {{-1.78585, 2.10707, 1.}, {10.4525, 8.22626, 1.}, {1., 0., 2.}}}

Created by Mathematica  (April 1, 2005) Valid XHTML 1.1!