(*^ ::[paletteColors = 128; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L3, e8, 24, "New York"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "New York"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 14, "New York"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 14, "New York"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 12, "New York"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 10, "New York"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "New York"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; ; fontset = name, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = Left Header, inactive, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; ; fontset = Left Footer, inactive, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Geneva"; ; fontset = clipboard, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = completions, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special1, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special2, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, center, M7, 12, "New York"; ; fontset = special3, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, right, M7, 12, "New York"; ; fontset = special4, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = special5, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ;] :[font = text; inactive; locked; preserveAspect; center; ] Numerical Methods: Mathematica Notebooks (c) John H. Mathews, 1996 To accompany the text: Numerical Methods: for Mathematics, Science, and Engineering, 2nd Ed, 1992 Prentice Hall, Inc. Englewood Cliffs, New Jersey 07632 U.S.A. Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6 Prentice Hall, International Editions: ISBN 0-13-625047-5 This free software is compliments of the author. E-mail address: in%"mathews@fullerton.edu" http://titan.fullerton.edu/~mathews/ ;[s] 6:0,1;67,0;92,1;157,3;159,1;232,2;486,-1; 4:1,17,12,New York,0,12,0,0,0;3,23,17,New York,0,18,0,0,0;1,19,14,New York,0,14,0,0,0;1,27,18,New York,32,14,0,0,0; :[font = text; inactive; locked; preserveAspect; center; ] CONTENTS ;[s] 2:0,1;8,0;9,-1; 2:1,17,12,New York,0,12,0,0,0;1,23,17,New York,1,18,0,0,0; :[font = text; inactive; locked; preserveAspect; ] Chapter 3. The Solution of Linear Systems AX = B Algorithm 3.1 Back Substitution Algorithm 3.2 Upper-Triangularization Followed by Back Substitution Algorithm 3.3 PA = LU Factorization with Pivoting Algorithm 3.4 Jacobi Iteration Algorithm 3.5 Gauss-Seidel Iteration ;[s] 2:0,1;49,0;278,-1; 2:1,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 3.1 (Back Substitution). To solve the upper-triangular system: a1,1 x1 + a1,2 x2 + ... + a1,n-1 xn-1 + a1,n xn = b1 a2,2 x2 + ... + a2,n-1 xn-1 + a2,n xn = b2 : : : an-1,n-1 xn-1 + an-1,n xn = bn-1 an,n xn = bn Proceed with the method only if all the diagonal elements are non-zero. First compute xn = bn/an,n and then use the rule xr = ÊÊ(br - ar,r+1 xr+1 - ... - a2,n-1 xn-1 - a2,n xn)/ar,rÊ for r = n-1,n-2, ...,1. Remark. Additional steps have been placed in the algorithm which perform the side calculation of finding the determinant of A. Section 3.3 Upper-Triangular Linear Systems Page 145 ;[s] 78:0,1;34,0;75,2;78,0;80,2;81,0;85,2;88,0;90,2;91,0;101,2;106,0;108,2;111,0;115,2;118,0;120,2;121,0;127,2;128,0;144,2;147,0;149,2;150,0;160,2;165,0;167,2;170,0;174,2;177,0;179,2;180,0;186,2;187,0;290,2;297,0;299,2;302,0;307,2;312,0;315,2;316,0;320,2;323,0;378,2;381,0;384,2;385,0;390,2;392,0;483,2;484,0;488,2;489,0;491,2;494,0;520,2;521,0;531,2;532,0;536,2;541,0;543,2;546,0;556,2;561,0;563,2;566,0;570,2;573,0;575,2;576,3;578,0;579,2;582,0;745,1;746,0;749,1;805,-1; 4:38,17,12,New York,0,12,0,0,0;3,17,12,New York,1,12,0,0,0;36,20,13,New York,64,9,0,0,0;1,25,16,New York,64,12,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[BackSub,a0,b0]; BackSub[a0_,b0_] := Module[{a,b,x,k,sum,mna,nb,ma,na}, a = a0; b = b0; mna = Dimensions[a]; nb = Length[b]; ma = mna[[1]]; na = mna[[2]]; det = 1; For[k=1,k<=ma,k++, det = det a[[k,k]]; ]; If[ma==na && det!=0, Module[{i,j}, x = Table[0,{i,ma}]; x[[ma]] = b[[ma]]/a[[ma,ma]]; Clear[i,j,k]; For[i=ma-1,1<=i,i--, sum = 0; For[j=i+1,j<=ma,j++, sum = sum + a[[i,j]] x[[j]] ]; x[[i]] = (b[[i]] - sum)/a[[i,i]]; ]; ], x = MatrixIsSingular]; Return[x]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 3, Example 3.12, Page 144. Use back substitution to solve the upper-triangular linear system AX = B: ;[s] 3:0,0;1,1;35,0;113,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] A = {{3,-1,2,3},{0,-2,7,-4},{0,0,6,5},{0,0,0,3}}; A = N[A,19]; TableForm[A] :[font = input; dontPreserveAspect; ] B = {20,-7,4,6}; B = N[B,19]; TableForm[B] :[font = input; dontPreserveAspect; ] X = BackSub[A,B]; TableForm[X] :[font = text; inactive; preserveAspect; ] Verify by matrix multiplication to see if AX - B = 0. :[font = input; preserveAspect; ] Y = A.X-B; TableForm[Y] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in routine. ;[s] 3:0,0;14,1;25,0;47,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; ] Z = LinearSolve[A,B]; TableForm[Z] TableForm[A.Z-B] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 3.2 (Upper-Triangularization Followed by Back Sub.) To construct the solution to AX = B, by first reducing the augmented matrix [A,B] to upper triangular form and then performing back substitution. Section 3.4 Gaussian Elimination and Pivoting Page 156 ;[s] 11:0,1;61,0;92,1;94,0;97,1;98,0;140,1;141,0;142,1;143,0;209,1;267,-1; 2:5,17,12,New York,0,12,0,0,0;6,17,12,New York,1,12,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] Let the coefficients of A and the constants of B be stored in the augmented matrix [A,B] = (ai,j) which has n rows and n+1 columns; that is, the column vector B = (ai,n+1) is stored in column n+1 of the augmented matrix [A,B]. The system can be written a1,1 x1 + a1,2 x2 +...+ a1,n-1 xn-1 + a1,n xn = a1,n+1 a2,1 x1 + a2,2 x2 +...+ a2,n-1 xn-1 + a2,n xn = a2,n+1 . . . . . : : : : : an-1,1 x1 + an-1,2 x2 +...+ an-1,n-1 xn-1 + an-1,n xn = an-1,n+1 an,1 x1 + an,2 x2 + ... + an,n-1 xn-1 + an,n xn = an,n+1 Row operations will be used to eliminate xp in column p. REMARKS. Rather than actually interchange rows, row interchanges are recorded in the pointer vector Row(p). If either the matrix A or vector B is to be used later, then a copy must be stored someplace else besides the augmented matrix. ;[s] 95:0,0;24,1;25,0;47,1;48,0;85,1;86,0;87,1;88,0;94,2;97,0;161,1;162,0;167,2;172,0;223,1;224,0;225,1;226,0;258,2;261,0;263,2;264,0;270,2;273,0;275,2;276,0;284,2;289,0;291,2;294,0;298,2;301,0;303,2;304,0;310,2;315,0;319,2;322,0;324,2;325,0;329,2;332,0;334,2;335,0;343,2;348,0;350,2;353,0;358,2;361,0;363,2;364,0;370,2;375,0;498,2;503,0;505,2;506,0;510,2;515,0;517,2;518,0;526,2;533,0;535,2;538,0;542,2;547,0;549,2;550,0;556,2;563,0;566,2;569,0;571,2;572,0;578,2;581,0;583,2;584,0;596,2;601,0;603,2;606,0;612,2;615,0;617,2;618,0;624,2;629,0;674,2;675,0;832,1;833,0;927,-1; 3:48,17,12,New York,0,12,0,0,0;8,17,12,New York,1,12,0,0,0;39,20,13,New York,64,9,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[GaussJordan,a0]; GaussJordan[a0_] := Module[{a,i,j,mna,n,n1,p,row,x, cond,Nonsingular,Singular}, a = a0; mna = Dimensions[a]; n = mna[[1]]; n1 = mna[[2]]; x = Table[0,{i,n}]; cond = Nonsingular; Clear[row]; For[j=1,j<=n,j++, row[j]=j]; For[p=1,p<=n-1,p++, Block[{k}, For[k=p+1,k<=n,k++, If[Abs[a[[row[k],p]]] > Abs[a[[row[p],p]]], t=row[p]; row[p]=row[k]; row[k]=t ]; ]; If[a[[row[p],p]] == 0, Print["A(",row[p],",",p,") = ",a[[row[p],p]]]; cond = Singular; Print["The matrix is singular."]; Goto[getout]; ]; Clear[k]; For[k=p+1,k<=n,k++, Module[{c}, m = a[[row[k],p]]/a[[row[p],p]]; For[c=p+1,c<=n+1,c++, a[[row[k],c]]=a[[row[k],c]]- m a[[row[p],c]] ]; ]; ]; ]; ]; If[a[[row[n],n]] == 0, Print["A(",row[n],",",n,") = ",a[[row[n],n]]]; cond = Singular; Print["The matrix is singular."]; Goto[getout]; ]; x[[n]] = a[[row[n],n+1]]/a[[row[n],n]]; Clear[k]; For[k=n-1,1<=k,k--, sum:= 0; For[c=k+1,c<=n,c++, sum = sum + a[[row[k],c]] x[[c]]; x[[k]] = (a[[row[k],n+1]]-sum)/a[[row[k],k]]; ]; ]; Label[getout]; Return[x]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 3, Exercise 8, Page 158. Use upper-triangularization followed by back-substitution to solve the linear system AX = B: ;[s] 3:0,0;1,1;33,0;130,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] A = {{4,8,4,0},{1,5,4,-3}, {1,4,7,2},{1,3,0,-2}}; A = N[A,19]; TableForm[A] :[font = input; dontPreserveAspect; ] B = {8,-4,10,-4}; B = N[B,19]; TableForm[B] :[font = input; preserveAspect; ] AB = {{4,8,4,0,8},{1,5,4,-3,-4}, {1,4,7,2,10},{1,3,0,-2,-4}}; AB = N[AB,19]; TableForm[AB] :[font = input; dontPreserveAspect; ] X = GaussJordan[AB]; TableForm[X] :[font = text; inactive; preserveAspect; ] Verify by matrix multiplication to see if AX - B = 0. :[font = input; dontPreserveAspect; ] Y = A.X - B; TableForm[Y] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in routine. ;[s] 3:0,0;14,1;25,0;47,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; ] Z = LinearSolve[A,B]; TableForm[Z] TableForm[A.Z-B] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 3.3 (PA = LU Factorization with Pivoting). To construct the solution to the linear system AX = B by performing the steps: 1. Find a permutation matrix P , lower-triangular matrix L and upper-triangular matrix U which satisfy PA = LU. 2. Compute PB and form the equivalent linear system LUX = PB. 3. Solve the lower-triangular system LY = PB for Y. 4. Solve the upper-triangular system UX = Y for X. Remarks: Since the diagonal elements of L are all 1's these values do not need to be stored. The coefficients of L below the main diagonal and the nonzero coefficients of U overwrite the matrix A. The permutation matrix P is not actually needed. The pointer vector row serves in its place to keep track of row interchanges. If matrix A is to be used later, then a copy must be stored someplace else. The value det is the determinant of A. Section 3.6 Triangular Factorizaton Page 175 ;[s] 51:0,1;54,0;101,1;103,0;106,1;107,0;164,1;165,0;193,1;194,0;230,1;231,0;249,1;251,0;254,1;256,0;270,1;272,0;314,1;317,0;320,1;322,0;361,1;363,0;366,1;368,0;375,1;376,0;415,1;417,0;420,1;421,0;428,1;429,0;471,1;472,0;545,1;546,0;603,1;604,0;626,1;627,0;652,1;653,0;698,1;701,0;768,1;769,0;874,1;875,0;877,1;925,-1; 2:25,17,12,New York,0,12,0,0,0;26,17,12,New York,1,12,0,0,0; :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] Off[General::spell1]; Clear[LUfactor,a0,row]; LUfactor[a0_] := Module[{a,i,j,p,mna,n,m,rowk,rowp}, a = a0; mna = Dimensions[a]; n = mna[[1]]; m = mna[[2]]; det = 1; For[j=1,j<=n,j++, row[j]=j]; For[p=1,p<=n-1,p++, Block[{c,k}, For[k=p+1,k<=n,k++, If[Abs[a[[row[k],p]]]>Abs[a[[row[p],p]]], t=row[p]; row[p]=row[k]; row[k]=t; det=-det]; ]; det = det a[[row[p],p]]; If[det == 0, Print["The matrix is singular."]; Goto[getout]; ]; Clear[k]; For[k=p+1,k<=n,k++, rowk = row[k]; rowp = row[p]; a[[rowk,p]] = a[[rowk,p]]/a[[rowp,p]]; For[c=p+1,c<=n,c++, a[[rowk,c]] = a[[rowk,c]] - a[[rowk,p]] a[[rowp,c]];];];];]; Label[getout]; det = det a[[row[n],n]]; Return[a]; ]; :[font = input; dontPreserveAspect; endGroup; ] Clear[SolveLU,lu0,b0]; SolveLU[lu0_,b0_] := Module[{a,c,i,j,k,n,p,x,rowk,sum}, a = lu0; b = b0; n = Length[a]; x = Table[0,{i,n}]; For[k=1,k<=n,k++, If[a[[row[k],k]]==0, Goto[getout] ]; ]; x[[1]] = b[[row[1]]]; Clear[k]; For[k=2,k<=n,k++, sum = 0; rowk = row[k]; For[c=1,c<=k-1,c++, sum = sum + a[[rowk,c]] x[[c]] ]; x[[k]] = b[[rowk]] - sum; ]; x[[n]] = x[[n]]/a[[row[n],n]]; Clear[k]; For[k=n-1,1<=k,k--, sum = 0; rowk = row[k]; For[c=k+1,c<=n,c++, sum = sum + a[[rowk,c]] x[[c]] ]; x[[k]] = (x[[k]] - sum)/a[[rowk,k]]; ]; Label[getout]; Return[x]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 3, Exercise 7, Page 178. Use LU factorization with pivoting to solve the linear system AX = B: ;[s] 3:0,0;1,1;33,0;107,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] A = {{2,4,-4,0}, {1,5,-5,-3}, {2,3,1,3}, {1,4,-2,2}}; A = N[A,19]; TableForm[A] :[font = input; dontPreserveAspect; ] LU = LUfactor[A]; TableForm[LU] :[font = input; dontPreserveAspect; ] B = {12,18,8,8}; B = N[B,19]; TableForm[B] :[font = input; dontPreserveAspect; ] X = SolveLU[LU,B]; TableForm[X] :[font = text; inactive; preserveAspect; ] Verify by matrix multiplication to see if AX - B = 0. :[font = input; dontPreserveAspect; ] Y = A.X - B; TableForm[Y] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in routine. ;[s] 3:0,0;14,1;25,0;47,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; ] Z = LinearSolve[A,B]; TableForm[Z] TableForm[A.Z-B] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 3.4 (Jacobi Iteration). To solve the linear system AX = B by starting with P0 = 0, and generating a sequence {Pk} which converges to the solution P, i.e. AP = B. A sufficient condition for the method to be applicable is that A is diagonally dominant. Section 3.7 Iterative Methods for Linear Systems Page 186 ;[s] 23:0,1;33,0;61,1;63,0;66,1;67,0;87,1;88,2;89,0;92,1;93,0;122,1;123,2;124,0;159,1;160,0;169,1;171,0;174,1;175,0;243,1;244,0;272,1;333,-1; 3:10,17,12,New York,0,12,0,0,0;11,17,12,New York,1,12,0,0,0;2,20,13,New York,64,9,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[Jacobi,A,B,V,P0,Tol,max]; Jacobi[A_,B_,P0_,Tol_,max_] := Block[{r,c,k,n,cond,stat,X,Xold,row, DiagonallyDominant,NotDominant}, norm[V_] := Sqrt[V^2/.List->Plus]//N; sep = 1; k = 0; cond = 1; stat = DiagonallyDominant; X = P0 //N; Xold = P0; n = Length[X]; For[r=1,r<=n,r++, row = 0; For[c=1,c<=n,c++, row = row + Abs[A[[r,c]]]//N]; If[row >= 2 Abs[A[[r,r]]], cond = 0] ]; If[cond == 0, stat=NotDominant; Goto[getout]]; While[k Tol, For[r=1,r<=n,r++, sum = B[[r]]//N; For[c=1,c<=n,c++, If[c!=r,sum = sum-A[[r,c]]Xold[[c]]//N];]; X[[r]] = sum/A[[r,r]]//N ]; sep = norm[X - Xold]; Print[N[X,14]]; Xold = X; k = k + 1]; Label[getout]; Return[N[X,14]]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 3, Example 3.25, Page 180. Use Jacobi iteration to solve the linear system AX = B: ;[s] 3:0,0;1,1;35,0;95,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] A = {{4,-1,1}, {4,-8,1}, {-2,1,5}}; A = N[A,19]; TableForm[A] :[font = input; dontPreserveAspect; ] B = {7,-21,15}; B = N[B,19]; TableForm[B] :[font = input; dontPreserveAspect; ] P0 = {1,2,2}; TableForm[P0] :[font = input; dontPreserveAspect; ] X = Jacobi[A,B,P0,0.00000001,100] TableForm[X] :[font = text; inactive; preserveAspect; ] Verify by matrix multiplication to see if AX - B = 0. :[font = input; preserveAspect; ] Y = A.X - B :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in routine. ;[s] 3:0,0;14,1;25,0;47,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; ] Z = LinearSolve[A,B]; TableForm[Z] TableForm[A.Z-B] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 3.5 (Gauss-Seidel Iteration). To solve the linear system AX = B by starting with P0 = 0, and generating a sequence {Pk} which converges to the solution P, i.e. AP = B. A sufficient condition for the method to be applicable is that A is diagonally dominant. Section 3.7 Iterative Methods for Linear Systems Page 187 ;[s] 23:0,1;39,0;67,1;69,0;72,1;73,0;91,1;92,2;93,0;96,1;97,0;127,1;128,2;129,0;165,1;166,0;175,1;177,0;181,1;182,0;250,1;251,0;279,1;340,-1; 3:10,17,12,New York,0,12,0,0,0;11,17,12,New York,1,12,0,0,0;2,20,13,New York,64,9,0,0,0; :[font = input; dontPreserveAspect; ] Off[General::spell1]; Clear[GaussSeidel,A,B,V,P0,Tol,max]; GaussSeidel[A_,B_,P0_,Tol_,max_] := Block[{r,c,k,n,X,Xold,cond,row,stat, DiagonallyDominant,NotDominant}, norm[V_] := Sqrt[V^2/.List->Plus]//N; sep = 1; k = 0; cond = 1; stat = DiagonallyDominant; X = P0 //N; Xold = P0; n = Length[X]; For[r=1,r<=n,r++, row = 0; For[c=1,c<=n,c++, row = row + Abs[A[[r,c]]]//N]; If[row >= 2 Abs[A[[r,r]]], cond = 0] ]; If[cond == 0, stat=NotDominant; Goto[getout]]; While[k Tol, For[r=1,r<=n,r++, sum = B[[r]]//N; For[c=1,c<=n,c++, If[c!=r, sum = sum-A[[r,c]]X[[c]]//N];]; X[[r]] = sum/A[[r,r]]//N ]; sep = norm[X - Xold]; Print[N[X,14]]; Xold = X; k = k + 1]; Label[getout]; Return[N[X,14]]; ]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 3, Example 3.25, Page 180. Use Gauss-Seidel iteration to solve the linear system: ;[s] 3:0,0;1,1;35,0;93,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = input; dontPreserveAspect; ] A = {{4,-1,1}, {4,-8,1}, {-2,1,5}}; A = N[A,19]; TableForm[A] :[font = input; dontPreserveAspect; ] B = {7,-21,15}; B = N[B,19]; TableForm[B] :[font = input; dontPreserveAspect; ] P0 = {1,2,2}; TableForm[P0] :[font = input; dontPreserveAspect; ] X = GaussSeidel[A,B,P0,0.00000001,100] TableForm[X] :[font = text; inactive; preserveAspect; ] Verify by matrix multiplication to see if AX - B = 0. :[font = input; preserveAspect; ] Y = A.X - B; TableForm[Y] :[font = text; inactive; preserveAspect; ] Compare with Mathematica's built in routine. ;[s] 3:0,0;14,1;25,0;47,-1; 2:2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0; :[font = input; dontPreserveAspect; ] Z = LinearSolve[A,B]; TableForm[Z] TableForm[A.Z-B] ^*)