Num. Math.—6

Iteráció

Mathematica-ban az n-edik iterált Nest[f,x0,n] 
Iteráltsorozat {x0, f[x0], ff[x0],...} NestList[f,x0,n]

Nest[f,x0,3]

Nest[f,x0,0]

Nest[f,x0,1]

NestList[f,x0,3]

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

Nest[f,1,3]

NestList[f,1,3]

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

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

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

Jacobi-Iteráció

Cél: megoldani az Ax=a lin. 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 eüh. mátrix ritka.

Szüks. és elegendõ felt. (Neumann-mátrix sor-konvergencia) spektrálsugár<1
Iteráció: Ax=a helyett x=Bx+b
Szüks. és elegendõ felt. (konvergencia) spektrálsugár<1

Elegendõ: valamelyik normára  ||B||<1

Jacobi

Ha a diagonális elemek nemzérók

A x=a⟹x=B x+b,ahol

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

Nemszinguláris együtthatómátrix.

Remove[A,B,a,b]

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

Det[A]

(Pontos) Megoldásvektor (Exact Solution)

ES1=LinearSolve[A,a]

Ellenôrzés

A.ES1

Iteráció elõkészítése

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!

(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

Feladat

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}}));

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].a;

f[x_]:=B.x+b

B

b=Flatten[b,1]

NestList[f,b,15]//N

Utolsó eleme a listának, utolsó iterált.

Last[%]

LinearSolve[A,a]//N

LinearSolve[A,a]//TableForm

GS-Iteráció

Feladat

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

Összevetés Jacobival.

Összevetés a Mathematica bép. megoldó rutinjával.

LinearSolve[A,a]//N

LinearSolve[A,a]

diag dom, partició

SDM[n_]:={PadRight[{4,-1},n],Sequence@@Table[ReplacePart[Table[0,{n}],{j-1->-1,j->4,j+1->-1}],{j,2,n-1}],PadLeft[{-1,4},n]}

Table[MatrixForm[SDM[n]],{n,3,6}]

Table[Det[SDM[n]],{n,3,10}]

Table[Eigenvalues[SDM[n]],{n,3,10}]

N[%]

PartJac[A_]:={UpperTriangularize[A,1],LowerTriangularize[A,-1],DiagonalMatrix[Table[A[[j,j]],{j,Length[A]}]]}

(A=Table[w[j,k],{j,4},{k,4}])//MatrixForm

MatrixForm[UpperTriangularize[A,4]]

Map[MatrixForm,PartJac[SDM[5]]]

PartJac[SDM[7]]

Map[MatrixForm,%]

Eigenvalues[-Inverse[%[[3]]].(%[[1]]+%[[2]])]

N[%]

GS expicit (n=3) [(Subscript[x, j])=((Subscript[x, j])^((m-1))),(Subscript[y, j])=((Subscript[x, j])^((m))) átnevezéssel]

-DiagonalMatrix[1/{a11,a22,a33}].({{0,a12,a13},{0,0,a23},{0,0,0}}.{y1,y2,y3}+{{0,0,0},{a21,0,0},{a31,a32,0}}.{x1,x2,x3}-{a14,a24,a34})

Jacobi explicit

%/.{y1->x1,y2->x2,y3->x3}

Példa ||A|| definíció alapján közelítve, ill mint felt. optimalizálási probléma megoldása

A=({{12,-2,1},{-2,15,3},{1,-2,10}});

T=Table[Table[Random[Real,{-1,1}],{3}],{10^5}];

Max[Map[Norm[A.#]/Norm[#]&,T]]

NMaximize[{Norm[A.{x1,x2,x3}],x1^2+x2^2+x3^2==1},{x1,x2,x3}]

Maximize[{Norm[A.{x1,x2,x3}],x1^2+x2^2+x3^2==1},{x1,x2,x3}]