Power Method

a = {{6, 5, -5}, {2, 6, -2}, {2, 
            5, -1}}; 
Eigensystem[a]

PowerMethod[a_, it_] := 
      Module[{lam, x, x0, y, phi}, n = Length[a[[1]]]; 
        x0 = Table[Random[], {i, n}];  
        phi[x_] := x[[3]]; 
Do[y = a . x0; 
          lam = phi[y]/phi[x0]; x = y/Max[Abs[y]]; Print[lam, "   ", x]; 
          x0 = x, {i, it}]; {lam, x}
]

pm = PowerMethod[a, 10];

pm[[1]]
pm[[2]]

Los izbor vektora x0

x0 = {1, 2, 3. };  phi[x_] := x[[3]];
Do[y = a . x0; lam = phi[y]/phi[x0]; x = y/Max[Abs[y]]; 
      Print[lam, "   ", x]; x0 = x, {i, 5}]

Zasto?

Eigensystem[a]

LinearSolve[{{1, 0, 1}, {1, 1, 0}, {1, 1, 1}}, {1, 2, 3}]

b = Inverse[a]
Eigensystem[b]

PowerMethod[b, 10];

(* Ako stavim phi[x]=x[[2]], PowerMethod za b ne konvergira prema najvecoj svojstvenoj vrijednosti. Zasto? *)

Inverse Power Method

x0 = {2. , 3, 4};
Do[y = LinearSolve[a, x0]; lam = phi[y]/phi[x0]; x = y/Max[Abs[y]]; 
      Print[lam, "   ", x]; x0 = x, {i, 10}]

Popravak

\(LU[a_, n_] := 
      Module[{i, j, k}, \ 
        L = Table[0, {i, n}, {j, n}]; \ \ U = 
          Table[0, {i, n}, {j, n}]; \ 
		Do[
          L[\([i, i]\)] = 1, {i, 
            n}]; 
	\ \ \ 	Do[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Do[
\ \ \ \ \ 	U[\([i, j]\)] = 
              a[\([i, j]\)] - 
                Sum[L[\([i, k]\)]\ U[\([k, j]\)], {k, i - 1}], 
				{j, 
              i, n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Do[
		L[\([j, 
                  i]\)] = \((a[\([j, i]\)] - 
                    Sum[L[\([j, k]\)] U[\([k, i]\)], {k, i - 1}])\)/
                U[\([i, i]\)], 
					{j, i + 1, n}]	
			, {i, 
            n}]; 
		Print["A=", MatrixForm[a], "  L=", 
          MatrixForm[L], ",  U=", MatrixForm[U]]; 
{L, 
          U}
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)
\(FS[n_, a_, b_] := 
      Module[{i, j}, \ \ \ x = Table[0, {i, n}]; 
		\ \ \ x[\([1]\)] = 
          b[\([1]\)]/
            a[\([1, 1]\)]; 
	\ \ \ 	Do[
		\ \ \ \ \ \ \ \ \ 	x[\([i]\)] = \((b[\([i]\)] - Sum[a[\([i, j]\)]\ x[\([j]\)], {j, i - 1}])\)/
              a[\([i, i]\)]
				\ \ \ \ \ \ \ , {i, 2, n}]; 
        x
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
 (*\ Algoritam\ BS\ *) \)
\(BS[n_, a_, b_] := 
      Module[{i, j, k}, \ \ \ \ \ \ \ x = 
          Table[0, {i, n}]; 
		x[\([n]\)] = 
          b[\([n]\)]/a[\([n, n]\)]; 
		Do[
          k = n - i; 
			\ \ \ \ \ x[\([k]\)] = \((b[\([k]\)] - 
                  Sum[a[\([k, j]\)]\ x[\([j]\)], {j, k + 1, n}])\)/
              a[\([k, k]\)]
\ \ \ \ 			\ , {i, 1, n - 1}]; 
        N[x, 5]
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

LU[a, n];

MODULI

\(\(\( (*\ Algoritam\ LU\ *) \)\(
\)\(LU[a_, n_] := 
      Module[{i, j, k}, \ 
        L = Table[0, {i, n}, {j, n}]; \ \ U = 
          Table[0, {i, n}, {j, n}]; \ 
		Do[
          L[\([i, i]\)] = 1, {i, 
            n}]; 
	\ \ \ 	Do[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Do[
\ \ \ \ \ 	U[\([i, j]\)] = 
              a[\([i, j]\)] - 
                Sum[L[\([i, k]\)]\ U[\([k, j]\)], {k, i - 1}], 
				{j, 
              i, n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Do[
		L[\([j, 
                  i]\)] = \((a[\([j, i]\)] - 
                    Sum[L[\([j, k]\)] U[\([k, i]\)], {k, i - 1}])\)/
                U[\([i, i]\)], 
					{j, i + 1, n}]	
			, {i, 
            n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {L, 
          U}
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]

    FS[n_, a_, b_] := 
      Module[{i, j}, \ \ \ x = Table[0, {i, n}]; 
		\ \ \ x[\([1]\)] = 
          b[\([1]\)]/
            a[\([1, 1]\)]; 
	\ \ \ 	Do[
		\ \ \ \ \ \ \ \ \ 	x[\([i]\)] = \((b[\([i]\)] - Sum[a[\([i, j]\)]\ x[\([j]\)], {j, i - 1}])\)/
              a[\([i, i]\)]
				\ \ \ \ \ \ \ , {i, 2, n}]; 
        x
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
 (*\ Algoritam\ BS\ *) 

    BS[n_, a_, b_] := 
      Module[{i, j, k}, \ \ \ \ \ \ \ x = 
          Table[0, {i, n}]; 
		x[\([n]\)] = 
          b[\([n]\)]/a[\([n, n]\)]; 
		Do[
          k = n - i; 
			\ \ \ \ \ x[\([k]\)] = \((b[\([k]\)] - 
                  Sum[a[\([k, j]\)]\ x[\([j]\)], {j, k + 1, n}])\)/
              a[\([k, k]\)]
\ \ \ \ 			\ , {i, 1, n - 1}]; 
        N[x, 5]
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
 (*\ 
      Power\ Method\ *) 

    PowerMethod[a_, it_] := 
      Module[{lam, x, x0, y, phi}, n = Length[a[\([1]\)]]; 
        x0 = Table[Random[], {i, n}]; \ 
        phi[x_] := x[\([3]\)]; 
Do[y = a . x0; 
          lam = phi[y]/phi[x0]; 
          x = y/Max[Abs[y]];  (*\ \(Print[lam, "   ", x];\)*) \ 
          x0 = x, {i, it}]; {lam, 
          x}
]

 (*\ 
      Inverse\ Power\ Method\ *) 

    InversePowerMethod[a_, it_] := 
      Module[{lam, x, x0, y, phi}, n = Length[a[\([1]\)]]; LU[a, n]; 
        x0 = Table[Random[], {i, n}]; \ 
        phi[x_] := x[\([3]\)]; 
Do[xx = FS[n, L, x0]; 
          y = BS[n, U, xx]; lam = phi[y]/phi[x0]; 
          x = y/Max[Abs[y]];  (*\ \(Print[lam, "   ", x];\)\ *) \ 
          x0 = x, {i, it}]; {lam, x}
]\)\)\)

Test primjer

n = 30 ;a = 
        Table[10 Random[], {i, n}, {j, n}] ;

pm = PowerMethod[a, 10];

pm[[1]]

ipm = InversePowerMethod[a, 10];

ipm[[1]]