Easy Sunday...Iterative Methods : Ax=b

[toshk]

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