Jeremy Knight Math 639 Algorithm 7.1: Jacobi Iterative
<Text-field style="Heading 1" layout="Heading 1">jacobi_iterate(n,A,b,x_0,tol,N,displa)</Text-field> jacobi_iterate := proc(n,A,b,x_0,tol,N,display) #n: scalar number of equations and unknown; #A: Coefficient array #b: Constant Vector #XO: Initial approximation vector #tol: tolerance #N: Max iterations local i, j, k, m, x,T,X0,x_k, DX, euclid_norm; #initialize a matrix of vectors to store the iterates #option trace; #X0:=Matrix(); #X0(1,1):=x_0; X0:=x_0; k:=1; while k <= N do #initialize a vector x to store the approximation values during the loop x:=Array(1..n); if display=1 then print(k-1,evalf(X0)); end if; x_k:=X0; for i from 1 to n do T:=0: for j from 1 to n do if i<>j then T:=T+(A(i,j)*x_k(j)): end if: end do; x(i):=(1/(A(i,i)))*((-1*T)+b(i)): end do: X0:=x: #Calculate the norm of ||x_(k+1)-x_k|| DX:=abs(X0 - x_k); euclid_norm:=DX(1); for m from 2 to n do if DX(m)>euclid_norm then euclid_norm:=DX(m); end if; end do; #Compare Euclidean Norm to tolerance if euclid_norm < tol then print("tolerance reached"); print(k,evalf(X0)); break; end if; if k>N then break; end if; k:=k+1; end do; if k=N+1 then print("tolerance not reached"); print(k-1,evalf(X0)); end if; end proc: