Easy Sunday...Iterative Methods : Ax=b - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: HP Prime (/forum-5.html) +--- Thread: Easy Sunday...Iterative Methods : Ax=b (/thread-9036.html) Easy Sunday...Iterative Methods : Ax=b - toshk - 09-10-2017 06:04 PM Stationary [[A]] Square [[A]] Can’t afford inversion of [[A]] [[A]] strictly diagonally dominant CAS programming: with Max Iterations of 10000. Jacobi Iteration: JbI([[A]],[guess],[b],error) Gauss-Seidel iteration: GsI([[A]],[guess],[b],error) Steepest Decent Method: SdI([[A]],[guess],[b],error) Conjugate Gradient Method: CgI([[A]],[guess],[b],error) example: JbI([[0.1,0,0],[0,0.1,−0.1],[0,−0.1,0.1]],[[0],[0],[0]],[[0.1],[−0.1],[0.1]],0.0001) Code: ``` #pragma mode( separator(.,;) integer(h32) ) #cas JbI(fun,k,b,ee):= BEGIN LOCAL r,d,T1,kk,d,dii; LOCAL roew,I; I:=0; T1:=1; kk:=[]; d:=diag(diag(fun)); dii:=diag(diag(fun) .^ -1); r:=fun-d; roew:=rowDim(fun); WHILE T1>ee DO kk:=(dii)*(b-(r*k)); IF Zr==1823 THEN kk:=dii* (-l*kk-u*k+b); END; I:=I+1; T1:=ABS(fun*kk-b);//residual error //IF ABS(k-kk)<=e THEN BREAK; END; //error IF I==10000 THEN BREAK; END;//forced out k:=kk; END; IF (T1>ee) AND (I==10000) THEN msgbox("Try Other Methods"); END; Zr:=0; return kk; END; #end #cas GsI(fun,k,b,ee):= BEGIN LOCAL Zr,UD,c,d,l,u; Zr:=1823; UD:=matrix(colDim(fun),colDim(fun)); c:=MAKELIST(X,X,1,colDim(fun)); FOR d FROM 1 TO length(c) DO pchD(c[d]); END; JbI(fun,k,b,ee); END; #end #cas pchD(dd) BEGIN LOCAL bb; FOR bb FROM dd TO colDim(fun) DO UD[dd,bb]:=fun[dd,bb]; END; l:=fun-UD; u:=UD-diag(diag(fun)); END; #end #cas SdI(fun,k,b,ee):= BEGIN LOCAL ro,z,al,I; I:=0; T1:=1; ro:=b-fun*k; WHILE ((ABS(ro)/ABS(b)) > ee) AND (I < 10000) DO z:=fun*ro; al:=(transpose(ro)*ro) ./ (transpose(ro)*z); k:=k+al[1,1]*ro; ro:=ro-al[1,1]*z; I:=I+1; END; return (k); END; #end #cas CgI(fun,k,b,ee):= BEGIN LOCAL ro,z,al,I; I:=0; po:=b-fun*k; ro:=po; WHILE ((ABS(ro)/ABS(b)) > ee) AND (I < 10000) DO z:=fun*po; al:=(transpose(po)*ro) ./ (transpose(z)*po); k:=k+al[1,1]*po; ro1:=ro-al[1,1]*z; bl:=(transpose(ro1)*ro1) ./ (transpose(ro)*ro); po:=ro1+bl[1,1]*po; I:=I+1; ro:=ro1; END; return (k); END; #end```