program endgueltig; { Attention: This program only runs for } { interval H-matrices } { with 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, sum1,sum2,klammer : 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; {calculate a crude enclosure of [x*]} {first 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]); {then, 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); { jetzt berechnen wir |[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; writeln; zaehl:=0; repeat zaehl:= zaehl+1; y_alt:=y_neu; for i:= 1 to n do begin if i=1 then begin sum1:=0; for k:=2 to n do sum1:=sum1 + M[1,k]*y_alt[k]; klammer:= -sum1 - q[1]; end else if i=n then begin sum2:=0; for k:=1 to n-1 do sum2:=sum2+M[n,k]*y_neu[k]; klammer:= -sum2-q[n]; end else begin sum1:=0; sum2:=0; for k:=i+1 to n do sum1:=sum1+M[i,k]*y_alt[k]; for k:= 1 to i-1 do sum2:= sum2+M[i,k]*y_neu[k]; klammer:=-sum1-sum2-q[i]; end; y_neu[i]:=imax(klammer/M[i,i])**y_alt[i]; end; 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.