{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "Highlight" -1 256 "" 0 0 0 255 0 1 0 1 0 0 0 0 
0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 
1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1
" -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }3 
1 0 0 8 4 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 256 1 {CSTYLE "" -1 
-1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 
0 1 }{PSTYLE "Heading 1" -1 257 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 
1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 }}
{SECT 0 {PARA 3 "" 0 "" {TEXT -1 55 "Solving Systems of Linear Equatio
ns:\nIterative Methods\n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "
restart;\nwith(linalg):" }}}{SECT 1 {PARA 256 "" 0 "" {TEXT -1 16 "A S
imple Example" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "A := matrix
([[2,1],[1,2]]);\nb := vector([6,6]);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 16 "Jacobi iteration" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 228 "x
 := vector([0,0]);\nfor k from 1 to 22 do \n   u1 := (b[1]-A[1,2]*x[2]
)/A[1,1];\n   u2 := (b[2]-A[2,1]*x[1])/A[2,2];\n   x[1] := evalf(u1):
\n   x[2] := evalf(u2):\n   printf(\"k = %2d:   x[1] = %f,   x[2] = %f
\\n\", k, x[1], x[2]); \nod:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "G
auss-Seidel iteration" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 202 "x := vect
or([0,0]);\nfor k from 1 to 12 do \n   x[1] := evalf((b[1]-A[1,2]*x[2]
)/A[1,1]);\n   x[2] := evalf((b[2]-A[2,1]*x[1])/A[2,2]);\n   printf(\"
k = %2d:   x[1] = %f,   x[2] = %f\\n\", k, x[1], x[2]); \nod:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 257 "" 0 
"" {TEXT -1 22 "The General Algorithms" }}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 84 "A symmetric positive, diagonally dominant matrix (all thr
ee methods should converge)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 175 "n :
= 10:\nA := 1/100*randmatrix(n,n):\nA := evalf(evalm(transpose(A)&*A +
 diag(seq(n, i=1..n)))):\nb := randvector(n):\nsol := evalf(linsolve(A
,b));\neigs := evalf(eigenvalues(A));" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 84 "A symmetric positive, but NOT diagonally dominant matrix \+
(Jacobi might not converge)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "n :
= 10:\nA := 1/10*randmatrix(n,n):\nA := evalm(transpose(A)&*A):\nb := \+
randvector(n):\nsol := evalf(linsolve(A,b));\neigs := evalf(eigenvalue
s(A));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 16 "Jacobi iteration" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 401 "Jacobi := proc(A, b, n, M, x, tol)
\n   local k, i, u, err;\n   u := vector(n, 0);\n   err := 99999;\n   \+
for k from 1 to M while err > tol do\n      for i from 1 to n do\n    \+
     u[i] := evalf((b[i] - sum(A[i,j]*x[j], j=1..i-1) - sum(A[i,j]*x[j
], j=i+1..n))/A[i,i]);\n      od;\n      for i from 1 to n do\n       \+
  x[i] := u[i];\n      od;\n      err := norm(x-sol);\n   od;\n   prin
t(k);\n   return evalm(x);\nend:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 76 "x := vector(n, 0);\nx := Jacobi(A, b, n, 1000, x, 0.0
01);\n'sol' = evalm(sol);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 22 "Gauss-Seidel iteration" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
321 "GaussSeidel := proc(A, b, n, M, x, tol)\n   local k, i, err;\n   \+
err := 99999;\n   for k from 1 to M while err > tol do\n      for i fr
om 1 to n do\n         x[i] := evalf((b[i] - sum(A[i,j]*x[j], j=1..i-1
) - sum(A[i,j]*x[j], j=i+1..n))/A[i,i]);\n      od;\n      err := norm
(x-sol);\n   od;\n   print(k);\n   return evalm(x);\nend:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "x := vector(n, 0);\nx := GaussSeide
l(A, b, n, 1000, x, 0.001);\n'sol' = evalm(sol);" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 31 "Gauss-Seidel iteration with SOR" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 428 "SOR := proc(A, b, n, M, x, omega, tol)\n   local \+
k, i, u, err;\n   u := vector(n,0);\n   err := 99999;\n   for k from 1
 to M while err > tol do\n      for i from 1 to n do\n         u[i] :=
 evalf( (1-omega)*x[i] + omega*(b[i] - sum(A[i,j]*u[j], j=1..i-1) - su
m(A[i,j]*x[j], j=i+1..n))/A[i,i]);\n      od;\n      for i from 1 to n
 do\n         x[i] := u[i];\n      od;\n      err := norm(x-sol);\n   \+
od;\n   print(k);\n   return evalm(x);\nend:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 79 "x := vector(n, 0);\nx := SOR(A, b, n, 1000, x, 1
.25, 0.001);\n'sol' = evalm(sol);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}}}{MARK "0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 
1 }{PAGENUMBERS 0 1 2 33 1 1 }