program endgueltig; { Attentiom: This program only runs,} { if [M] is an interval H-matrix } { with positive diagonal entries } use i_ari, mvi_ari, mv_ari; function iga(A:imatrix;b:ivector):ivector[lb(b)..ub(b)]; var k,i,j : integer; x,y : ivector[lb(b)..ub(b)]; hilf : interval; begin for k:=lb(b) to ub(b)-1 do begin for i:= k+1 to ub(b) do begin for j:= k+1 to ub(b) do begin a[i,j]:=a[i,j] - a[k,j] * a[i,k] /a[k,k]; end; {jetzt wird die untere Dreiecksmatrix besetzt} a[i,k]:=a[i,k]/a[k,k]; end; end; {Loese nun Ly=b} for i:= lb(b) to ub(b) do begin hilf:=0; for j:= lb(b) to i-1 do begin hilf:= hilf + a[i,j]*y[j]; end; y[i]:=b[i]- hilf; end; {Loese nun Rx=y} for i:= ub(b) downto lb(b) do begin hilf:=0; for j:= i+1 to ub(b) do begin hilf:= hilf + a[i,j]*x[j]; end; x[i]:= (y[i] - hilf)/a[i,i]; end; for i:= lb(b) to ub(b) do iga[i]:=x[i]; end; function imax(a : interval):interval; begin if (a.sup < 0) or (a.sup=0) then imax:=0 else if (a.inf > 0) or (a.inf=0) then imax:=a else imax:=intval(0,a.sup); end; function betrag(a:interval):real; begin if a.inf > 0 then betrag:=a.sup else if a.sup < 0 then betrag:=-a.inf else if -a.inf < a.sup then betrag:=a.sup else betrag:=-a.inf; end; procedure haupt(n:integer;stop:integer); var i,j,k,h,l,zaehl : integer; Hilfmatrix : imatrix[1..n,1..n]; ep : rvector[1..n]; x_1,v,u,alpha : ivector[1..n]; M,R : imatrix[1..n,1..n]; q : ivector[1..n]; hvector : ivector[1..n]; y_alt,y_neu,ystart : ivector[1..n]; y : ivector[1..n]; Loese : ivector[1..2*n]; unten,oben : real; hilf : interval; begin writeln('Give the interval matrix [M] row-wise'); writeln('[M] has to be an H-matrix with positive diagonal entries'); for i:= 1 to n do for j:= 1 to n do read(M[i,j]); writeln('Give the interval vector [q]'); for i:= 1 to n do read(q[i]); writeln; R:= -M; for i:= 1 to n do R[i,i]:=0; {first we calculate a crude enclosure of [x*]} {calculate P(I-P)-1*alpha} {with alpha:=q([x_0],[x_1])} {it is P(I-P)-1*alpha = <[D]>-1*|[R]|*-1*<[D]>*alpha} {calculate [x_1], if [x_0]=0} for i:= 1 to n do x_1[i]:=imax(-q[i]/M[i,i]); {set alpha := q([x_0],[x_1])=q(0,[x_1])} {and multiply with <[D]>} for j:= 1 to n do alpha[j]:=x_1[j].sup*M[j,j].inf; { calculate <[M]>-1* alpha} for i := 1 to n do for j:= 1 to n do if i=j then Hilfmatrix[i,j]:=M[i,j].inf else Hilfmatrix[i,j]:=-betrag(M[i,j]); u:=iga(Hilfmatrix,alpha); { calculate |[R]|* <[M]>-1* alpha} for i := 1 to n do for j:= 1 to n do Hilfmatrix[i,j]:= betrag(R[i,j]); hvector:=Hilfmatrix*u; for i:= 1 to n do v[i]:=hvector[i]/M[i,i].inf; for i:= 1 to n do ep[i]:= v[i].sup; for i:= 1 to n do ystart[i] := intval(x_1[i].inf -< ep[i], x_1[i].sup +> ep[i] ); for i:= 1 to n do if ystart[i].inf<0 then ystart[i].inf:=0; y_neu:=ystart; for j:= 1 to n do begin unten :=y_neu[j].inf; oben :=y_neu[j].sup; writeln('[',unten:23:0:-1,',',oben:23:0:1,']'); end; writeln; zaehl:=0; repeat zaehl:= zaehl+1; y_alt := y_neu; hvector:= R*y_alt-q; for i:=1 to n do y_neu[i]:= imax(hvector[i]/M[i,i])**y_alt[i]; until ((zaehl=stop) or (y_neu=y_alt)); writeln; write('We had ',zaehl,' iteration steps'); writeln; writeln('Enclosure of LB : '); writeln; for j:= 1 to n do begin unten :=y_neu[j].inf; oben :=y_neu[j].sup; writeln('[',unten:23:0:-1,',',oben:23:0:1,']'); end; end; var n,stop : integer; begin writeln('Give the dimension of the LCP'); read(n); writeln('How many iteration steps should be done at most ?'); read(stop); haupt(n,stop); end.