Jacobi és Gauss-Seidel iteráció

Iteráció Mathematica-val

Nest, NestList, NestWhile
A Mathematicában elegánsan képezhetjük az iteráltakat.

Nest[F,x,3]

NestList[F,x,3]

f[x_]:=x^2+1;

NestList[f,1,3]

Nest[f,1,3]

g[{x_,y_}]:={y,x+y}

NestList[g,{0,1},10]

NestList[g,{0,1},10][[All,1]]//TableForm

Spektrálsugár, mátrixnorma

Vektornormák

v={1,2,3};

Norm[v]

Norm[v,2]

Norm[v,1]

Norm[v,∞]

Mátrixnormák.
Itt: vektornormából származtatott mátrixnorma.
A default itt is a kettes.

Memo: ||A(Subscript[||, ∞])=(max_i)(∑_j)|(Subscript[a, i j])|
             ||A(Subscript[||, 1])=(max_j)(∑_i)|(Subscript[a, i j])|

A=(({{1, 9-3I, 4}, {2, 5, 6}, {-I, 6, 8}}));

Norm[A,1]

N[%]

Max[Map[Norm[A.#,1]/Norm[#,1]&,Table[Table[Random[Real,5],{3}],{10^6}]]]

Norm[A,∞]

B=(({{1, 2, 4}, {2, 5, 6}, {4, 6, 8}}));

{Norm[B,1],Norm[B,∞]}

A két norma (sor és oszlop) ugyanaz, mert a mátrix szimmetrikus

Kettes norma: spektrálnorma

Norm[B,2]

FullSimplify[%]

Max[Abs[Eigenvalues[B]]]

Ez nem teljesen véletlen: szimmetrikus mátrixok esetében a spektrálsugár megegyezik a 2-es normával, hiszen általában
||A(Subscript[||, 2])=(Sqrt[ρ((A^H)A)])

Ellenőrizzük ezt az A mátrixra

Norm[A,2]//N

ConjugateTranspose[A].A//MatrixForm

Sqrt[Max[Abs[Eigenvalues[ConjugateTranspose[A].A]]]//N]

minden normára ρ(A)≤||A||

Norm[A,1]//N

Norm[A,∞]

Abs[Eigenvalues[A]]//N

Neumann-sor

Konvergens-e? Ha igen, hova konvergál?

I+B+(B^2)+…+(B^m)+…=(∑_(m=0)^∞)(B^m)

B=(({{(1/2), 0}, {0, -(3/4)}}));
I2=IdentityMatrix[2];

SB={I2+B,I2+B+B.B,I2+B+MatrixPower[B,2],MatrixPower[B,3]};Map[MatrixForm[#]&,SB]

Map[MatrixForm[N[#]]&,SB]

MySB[M_,n_]:=(∑_(i=0)^(n-1))MatrixPower[M,i]

Table[MySB[B,n]//MatrixForm,{n,20}]//N

Inverse[IdentityMatrix[2]-B]//N//MatrixForm

B2=(({{(1/4), (1/2)}, {(1/2), -(1/4)}})); ∑(((B2))^n)=?

B2=(({{(1/4), (1/2)}, {(1/2), -(1/4)}}));

B2

Table[MySB[B2,n]//MatrixForm,{n,20,23}]//N

Norm[B2,1]

Inverse[IdentityMatrix[2]-B2]//N//MatrixForm

Inverse[IdentityMatrix[2]-B2]//MatrixForm

Cél: megoldani az Ax=a lineáris egyenletrendszert numerikusan. A megoldást fokozatosan közelítjük. Konvergencia kérdése, hatékonyság.
Tipikus alkalmazási terület: Parciális differenciálegyenletek numerikus megoldása. Ezeknél az együtthatómátrix ritka.

Jacobi-módszer: Ha a diagonális elemek nemzérók

A x=a⟹x=B x+b, ahol A=(L+D+U),

B=-(D^(-1))(L+U), b=(D^(-1))a

Ax=a x megoldásvektorának közelítése Jacobi-módszerrel, ahol:

A=(({{4, -1, 0}, {-1, 4, -1}, {0, -1, 4}}))a=(({{2}, {6}, {2}}))

Oldd meg a Jacobi iterációs módszerrel!

Nemszinguláris együtthatómátrix.

Clear[A,B,a,b]

A=(({{4, -1, 0}, {-1, 4, -1}, {0, -1, 4}}));a={2,6,2};

Det[A]

(Pontos) Megoldásvektor

ES1=LinearSolve[A,a]

Ellenőrzés:

A.ES1

Iteráció elõkészítés:

L={{0,0,0},{-1,0,0},{0,-1,0}};Di=DiagonalMatrix[{4,4,4}];
U={{0,-1,0},{0,0,-1},{0,0,0}};

Map[MatrixForm[#]&,{A,L,Di,U}]

MatrixForm[L+Di+U]

B=-Inverse[Di].(L+U);MatrixForm[B]

b=Inverse[Di].a

Eigenvalues[B]

Max[Abs[%//N]]

Konvergencia garantált!

Iterációs mag:

b

f1[x_]:=B.x+b

B, b mindig globális!

Iterációt megvalósító (rekurzív) algoritmus, Addig ismétli a szorzást és a vektor összeadást, amíg a második pm. le nem csökken 1-re.

(Subscript[x, 0])kezdõvektor tetszõlegesen választható, lehet pl. ugyanaz, mint b1, de lehet akár (Subscript[e, 1])={1,0,0} is.

x0={1,0,0};

B.x0+b

b

B//MatrixForm

f1[x0]

Nest[f1,x0,0]

Nest[f1,x0,1]

Nest[f1,x0,2]

Az ötödik iterált

Nest[f1,x0,5]//N

Az elsõ nyolc iterált

(T=NestList[f1,x0,8]//N)//ColumnForm

Hiba az n-edik iteráció után:

Map[Norm[#-ES1,1]&,T]//ColumnForm

A hibavektor normája gyorsan tart 0-hoz.

Hibabecslés:

||(Subscript[e, m])||≤||B||/(1-||B||) ||(Subscript[x, m])-(Subscript[x, m-1])|| (m=1,2,...)

Table[Norm[B,1]/(1-Norm[B,1])Norm[T[[j+1]]-T[[j]],1],{j,Length[T]-1}]//TableForm

Ax=a, ahol  A=(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))a=(({{3}, {-2}, {5}, {4}}))

Oldd meg a Jacobi iterációs módszerrel!

A=(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}));a=(({{3}, {-2}, {5}, {4}}));

Gauss-Seidel iteráció

A x=a⟹x=B x+b, ahol A=(L+D+U),

B=-(((D+L))^(-1))U, b=(((D+L))^(-1))a

Ax=a, ahol  A=(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))a=(({{3}, {-2}, {5}, {4}}))

Oldd meg a Gauss-Seidel iterációs módszerrel!

A=(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}));a=(({{3}, {-2}, {5}, {4}}));

Det[A]

Clear[Di,L,U]

L=(({{0, 0, 0, 0}, {2, 0, 0, 0}, {-1, 0, 0, 0}, {0, 2, -1, 0}}));

U=(({{0, -2, 1, 2}, {0, 0, 3, 1}, {0, 0, 0, 2}, {0, 0, 0, 0}}));

Di=DiagonalMatrix[{7,8,5,4}]

{L+U+Di//MatrixForm,A//MatrixForm}

B=-Inverse[Di+L].(U);b=Inverse[Di+L].a;

f[x_]:=B.x+b

B

b=Flatten[b,1]

NestList[f,b,8]//N

LinearSolve[A,a]//N

LinearSolve[A,a]