program Endversion;

{ f(x,s,t) = s + Arctan(xy') + pi }

use i_ari, mvi_ari, mv_ari;


{ Neue Idee:  waehle n=9 dann ist h=a/10  }
{ waehle dann N(k)=99 + k*100             }
{ dann gilt mit H:= a/((k+1)*100):        }
{ z[i]=i*H und x[j]= j*h                  }
{ wegen h=(k+1)*10*H gilt dann            }
{ x[j]=j*h=j*(k+1)*10*H                   }
{     =z[j*10*(k+1)]                      }



type  trigestalt = record 
                   diag,odiag,udiag : interval;
                   end;
      trimatrix = dynamic array[*] of trigestalt;


function f_kappa(y0:interval):interval;
var pihalbe:interval;
begin
pihalbe:= 2*Arctan(1);
f_kappa:=pihalbe;
end;

function f_Konst(y0:interval):interval;
var pi:interval;
begin 
pi:=4*Arctan(1);
f_Konst:=pi + y0; 
end;

function f_Lipp(y0:interval;a:real):interval;
begin
f_Lipp:=a;
end;


function f(x,s,t:interval):interval;
var pi :interval;
begin
pi:=4*Arctan(1);
f:=s+Arctan(x*t)+pi;
end; 

function f_x(x,s,t:interval): interval;
begin
f_x:=t/(1 + sqr(x*t));
end;

function f_s(x,s,t:interval):interval;
begin
f_s:=1;
end;


function f_t(x,s,t:interval):interval;
begin
f_t:=x/(1 + sqr(x*t));
end;






function G(x:rvector;ix:ivector;q:ivector):trimatrix[lbound(q)..ubound(q)];
var
argminy,argmaxz            : rvector[lbound(q)..ubound(q)];
i,j                        : integer;
nenner,zaehler,links       : interval;
wesen, hilf                : interval;


begin{function}


for i:= lbound(q) to ubound(q) do
begin{i=1(1)n}

{ 1. Schritt: Bestimme argminy }

for j := lbound(q) to ubound(q) do argminy[j]:=ix[j].sup;


if i=lbound(q) then hilf:=q[lbound(q)] - argminy[lbound(q)+1] * 0.5  
else
if i=ubound(q) then hilf:=q[ubound(q)] - argminy[ubound(q)-1] * 0.5 
else
hilf:=q[i] - argminy[i-1] * 0.5 - argminy[i+1] *0.5 ;

if not (hilf.inf < 0) then begin {Einheitsvektor uebergeben}
                         G[i].diag:=1;
                         if i = lb(q) then begin
                                           G[i].odiag:=0;
                                           G[i].udiag:=0;{kuenstlich}
                                           end
                         else if i = ub(q) then begin
                                                G[i].udiag:=0;
                                                G[i].odiag:=0;{kuenstlich}
                                                end
                              else begin
                                     G[i].odiag:=0;
                                     G[i].udiag:=0;
                                   end;
                       end
else 
begin{hilf.inf < 0 => 2. Schritt}

{ 2. Schritt: Bestimme argmaxz }

for j := lbound(q) to ubound(q) do argmaxz[j]:=ix[j].inf;


if i=lbound(q) then hilf:=q[lbound(q)] - argmaxz[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - argmaxz[ubound(q)-1] * 0.5 
else
hilf:=q[i] - argmaxz[i-1] * 0.5  - argmaxz[i+1] * 0.5 ;
if not (hilf.sup > 0) then begin {Zeile von M uebergeben}
                         G[i].diag:=1;
                         if i=lb(q) then begin
                                         G[i].odiag:=-0.5;
                                         G[i].udiag:=0;{kuenstlich}
                                         end
                         else if i=ub(q) then begin
                                              G[i].udiag:=-0.5;
                                              G[i].odiag:=0;{kuenstlich}
                                              end
                              else begin
                                     G[i].odiag:=-0.5;
                                     G[i].udiag:=-0.5;
                                   end;
                       end{2. Schritt}
else
begin{hilf.sup > 0 => 3. Schritt}

{ 3. Schritt }
if i=lbound(q) then hilf:=q[lbound(q)] - x[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - x[ubound(q)-1] * 0.5
else
hilf:=q[i] - x[i-1] * 0.5 - x[i+1] * 0.5 ;

if hilf.inf>0 then begin { 3.1}
               zaehler:=hilf;
               G[i].diag:=1;

               if i=lbound(q) then 
               begin
               hilf:= ( argminy[lbound(q)+1] - x[lbound(q)+1] ) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].odiag:=-0.5+wesen*0.5;
               G[i].udiag:=0;{kuenstlich}
               end
               else
               if i=ubound(q) then 
               begin
               hilf:= ( argminy[ubound(q)-1] - x[ubound(q)-1] ) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].udiag:=-0.5 + wesen*0.5;
               G[i].odiag:=0;{kuenstlich}
               end
               else
               begin
               hilf:=(argminy[i+1] - x[i+1] + argminy[i-1] - x[i-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].odiag:=-0.5 + wesen*0.5;
               G[i].udiag:=-0.5 + wesen*0.5;
               end;

               end {3.1}
else 
begin {3.2anfang}
if i=lbound(q) then hilf:=q[lbound(q)] - x[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - x[ubound(q)-1] * 0.5  
else
hilf:=q[i] - x[i-1] * 0.5 - x[i+1] * 0.5 ;
if hilf.sup < 0 then begin {3.2eigentlich}
               zaehler:=hilf;
               G[i].diag:=1;

               if i=lbound(q) then 
               begin
               hilf:=(argmaxz[lbound(q)+1] - x[lbound(q)+1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].odiag:=-wesen*0.5;
               G[i].udiag:=0;{kuenstlich}
               end
               else
               if i=ubound(q) then 
               begin
               hilf:=(argmaxz[ubound(q)-1] - x[ubound(q)-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].udiag:= -wesen*0.5;
               G[i].odiag:=0;{kuenstlich}
               end
               else
               begin
               hilf:=(argmaxz[i+1] - x[i+1] + argmaxz[i-1] - x[i-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               G[i].odiag:=-wesen*0.5;
               G[i].udiag:=-wesen*0.5;
               end;

     
               end {3.2eigentlich}
else
begin {3.3}
              wesen:=intval(0,1);
              G[i].diag:=1;
              if i=lbound(q) then begin
                                    G[i].odiag:= -wesen*0.5;
                                    G[i].udiag:=0;{kuenstlich}
                                  end
              else if i=ubound(q) then begin
                                         G[i].udiag:= -wesen*0.5;
                                         G[i].odiag:=0;{kuenstlich}
                                       end
                   else
                   begin
                    G[i].odiag:=-wesen*0.5;
                    G[i].udiag:=-wesen*0.5;
                   end;

end;{3.3}



end;{3.2anfang}

end;{3.Schritt}
end;{2.Schritt}

end;{i=1(1)n}
end;{function}

function igatri(A:trimatrix;b:ivector):ivector[lbound(b)..ubound(b)];

{ ------------------------------------------------------------------------- }

{ in diesem Programm kommen nur Tridiagonalmatrizen vor }

{ die LU Zerlegung reduziert sich auf }
{ for k=1 to n-1 do}
{ begin l[k+1,k]:=a[k+1,k]/a[k,k]; }
{       a[k+1,k+1]:= a[k+1,k+1] - l[k+1,k]*a[k,k+1]; }
{ end }

{ die Vorwaertszerlegung reduziert sich auf }
{ y[1]: b[1]; }
{ for i:= 2 to n do }
{ y[i]:=b[i] - l[i,i-1]*y[i-1] }

{ die Rueckwaertssubstitution reduziert sich auf }
{ x[n]:= y[n]/a[n,n]; }
{ for i:= n-1 downto 1 do }
{ x[i]:=(y[i]-a[i,i+1]*x[i+1])/a[i,i] }

{ man braucht also nur zwei Vektoren ltri und utri }

{ ------------------------------------------------------------------------- }

var ltri,utri,y,x : ivector[lbound(b)..ubound(b)];
    i,j,k         : integer;

begin
utri[lbound(b)]:= A[lb(b)].diag;
for k:= lbound(b) to ubound(b)-1 do 
         begin
              ltri[k+1]:= A[k+1].udiag/utri[k];
              utri[k+1]:= A[k+1].diag - ltri[k+1]*A[k].odiag;
         end;



    {Loese Ly=b}
     y[lbound(b)]:=b[lbound(b)];
     for i := 2 to ubound(b) do y[i]:=b[i]-ltri[i]*y[i-1];


     {Loese nun [R]x=y}

     x[ubound(b)]:= y[ubound(b)]/utri[ubound(b)];
     for i:= ubound(b)-1  downto lbound(b) do x[i]:= (y[i]-A[i].odiag*x[i+1])/utri[i];

for i:= lbound(b) to ubound(b) do igatri[i]:=x[i];

end;
   
function Hmin(x:rvector;q:ivector):ivector[lbound(q)..ubound(q)];
var 
    i            : integer;
    hilf         : interval;

begin

for i:= lbound(q) to ubound(q) do
begin

if i=lbound(q) then hilf:= x[lbound(q)] -  x[lbound(q)+1]*0.5 + q[lbound(q)]
else
   if i=ubound(q) then hilf:= x[ubound(q)] -  x[ubound(q)-1]*0.5 + q[ubound(q)]
   else hilf:= x[i] - x[i-1]*0.5 - x[i+1]*0.5 + q[i];
               
            
if x[i] > hilf.sup then Hmin[i]:=hilf
else
   if x[i] < hilf.inf then Hmin[i]:=x[i]
   else Hmin[i]:=intval(hilf.inf,x[i]);


end;
end;

function durchmesser(v : ivector):real;
var i            : integer;
    bisher, hilf : real;
begin
bisher:= diam(v[lbound(v)]);
for i:= lbound(v) + 1 to ubound(v) do
begin
     hilf:=diam(v[i]);
     if hilf > bisher then bisher:=hilf;
end;
durchmesser:=bisher;
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 drin(a,b : interval):boolean;
begin
if ((a.inf < b.inf) or (a.sup > b.sup))  then drin:=false
else drin:=true;
end;



function max( a,b : real):real;
begin
if a<b then max:= b
else        max:=a
end;

function betrag(a:real):real;
begin 
if a> 0 then betrag:= a
else betrag :=-a
end;


function leer(a,b : interval):boolean;
begin
if (a.inf <= b.sup) and (b.inf <= a.sup) then leer := false
else leer:= true;
end;


procedure mix(N,k:integer;y0:interval;var left:real;var a:real);
var

      
      qsesv, xstart, xneu, xalt : ivector[1..9];
      sum1,sum2,sum3,xhalb      : ivector[1..9];
      q,zstart,zalt,zneu        : ivector[1..N];
      schritt                   : interval;
      i,j, bis                  : integer;
      hilf,h1,h2                : interval;
      c,unten,oben,Fstrich      : real;
      arg1,arg2,arg3,part       : interval;
      Fausw,fx,fs,ft            : interval;
      kappa,Konst,Lipp          : interval;
      slopes                    : trimatrix[1..N];
      mitte                     : rvector[1..N];
      newton, Hminvector        : ivector[1..N];
      rechteseite               : ivector[1..N];
      count, zaehler, zaehl     : integer;
    
begin

kappa:=f_kappa(y0);
Konst:=f_Konst(y0);
Lipp:=f_Lipp(y0,a);

xstart:=intval(0,y0.sup);

schritt:=a/intval(10);

     arg1:=intval(0,a);
     arg2:=intval(0,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3);

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= 1 to 9 do qsesv[i]:=h1+h2;
     qsesv[1]:=qsesv[1]-y0;


zaehl:=0;
xneu:=xstart;

repeat  
     zaehl:= zaehl+1;


     xalt := xneu;
     

for i:= 1 to 9 do
   begin
  
      if i=1 then begin
	 sum2[1]:=0;
	 if zaehl > 1 then sum1[1]:=sum3[1]
         else sum1[1]:=xalt[2];
      end
      else if i=9 then begin
	 sum2[9]:=xhalb[8];
	 sum1[9]:=0;
      end
      else begin
	 if zaehl > 1 then sum1[i]:=sum3[i]
         else sum1[i]:=xalt[i+1];
         sum2[i]:=xhalb[i-1];
      end;
      xhalb[i]:=imax((sum1[i] + sum2[i] -qsesv[i])/2)**xalt[i];
   end;

   {jetzt xneu berechnen und verwende sum2[1..9]}

   for i:= 9 downto 1 do
   begin
      if i=9 then sum3[i]:=0
      else  begin
	 sum3[i]:=xneu[i+1];
      end;
     

      xneu[i]:=imax((sum2[i] + sum3[i] -qsesv[i])/2)**xhalb[i];

      
   end;



     if zaehl mod 10 = 0 then {neues [q]}
     begin
   
         bis:= 0;
         while (bis < 9) and (xneu[bis+1].inf > 0) do bis:= bis + 1;


     if bis > 1 then 
     begin
      
     for i:= 1 to bis-1 do
     begin
     arg1:=i*schritt;
     arg2:=xneu[i];
     hilf:=Lipp*arg2 + Konst*(a-arg1);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     if i> 1 then
     begin
     arg1:=intval((i-1)*<schritt.inf,(i+1)*>schritt.sup);
     arg2:=intval(xneu[i+1].inf,xneu[i-1].sup);
     hilf:=Lipp*xneu[i-1] + Konst*(a-(i-1)*schritt);
     arg3:=intval(-hilf.sup,0);
     end
     else
     begin
     arg1:=intval(0,2*>schritt.sup);
     arg2:=intval(xneu[2].inf,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);
     end;

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);
     part:=f(arg1,arg2,arg3);

     hilf:= abs(fx + fs * arg3 + ft * part );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/3;
     qsesv[i]:=h1+h2;
     end;
     qsesv[1]:=qsesv[1]-y0;



     arg1:=intval((bis-1)*<schritt.inf,a);
     arg2:=intval(0,xneu[bis-1].sup);
     hilf:=Lipp*xneu[bis-1] + Konst*(a-(bis-1)*schritt);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= bis to 9 do qsesv[i]:=h1+h2;
      
     end
     else writeln('no improvement possible yet');


     end;

 until xneu=xalt;


writeln('We have done ',zaehl,' iteration steps');








{ jetzt bestimmen wir zstart mit Hilfe von xneu }

for i:= 1 to 9 do
for j:= 1 to (k+1)*10 do
if i=1 then zstart[j]:= intval(xneu[1].inf,y0.sup)
else zstart[(i-1)*(k+1)*10+j]:=intval(xneu[i].inf,xneu[i-1].sup);

for j:=(k+1)*90+1 to N do
zstart[j]:=intval(0,xneu[9].sup);



zneu:=zstart;

     schritt:=a/intval(N+1);

     arg1:=intval(0,a);
     arg2:=intval(0,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3);

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;

     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= 1 to N do q[i]:=(h1+h2)/2;
     q[1]:=q[1]-y0/2;


zaehler:=0;

repeat
  zaehler:= zaehler+1;
  zalt:=zneu;
     
  for i:= 1 to N do mitte[i]:=zalt[i].inf;
       
  slopes:=G(mitte,zalt,q);

     count:=0;
     for i:= 1 to N do
     begin
          if not (diam(slopes[i].udiag)=0) then count:=count+1;
     end;
     writeln('slopes has ',count,' nondegenerate intervals');
 
  Hminvector:=Hmin(mitte,q);

  rechteseite:=igatri(slopes,Hminvector);
     
  newton:= mitte - rechteseite;

  for i:= 1 to N do zneu[i]:= zalt[i]**newton[i];


     if zaehler mod 10 = 0 then {neues [q]}
     begin
   
         bis:= 0;
         while (bis < N) and (zneu[bis+1].inf > 0) do bis:= bis + 1;


     if bis > 1 then 
     begin
      
     for i:= 1 to bis-1 do
     begin
     arg1:=i*schritt;
     arg2:=zneu[i];
     hilf:=Lipp*arg2 + Konst*(a-arg1);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     if i> 1 then
     begin
     arg1:=intval((i-1)*<schritt.inf,(i+1)*>schritt.sup);
     arg2:=intval(zneu[i+1].inf,zneu[i-1].sup);
     hilf:=Lipp*zneu[i-1] + Konst*(a-(i-1)*schritt);
     arg3:=intval(-hilf.sup,0);
     end
     else
     begin
     arg1:=intval(0,2*>schritt.sup);
     arg2:=intval(zneu[2].inf,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);
     end;

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);
     part:=f(arg1,arg2,arg3);

     hilf:= abs(fx + fs * arg3 + ft * part );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/3;
     q[i]:=(h1+h2)/2;
     end;
     q[1]:=q[1]-y0/2;

     arg1:=intval((bis-1)*<schritt.inf,a);
     arg2:=intval(0,zneu[bis-1].sup);
     hilf:=Lipp*zneu[bis-1] + Konst*(a-(bis-1)*schritt);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= bis to N do q[i]:=(h1+h2)/2;
     
     end
     else writeln('no improvement possible yet');


     end;

         
    
until zalt=zneu;


for i:= 1 to N do
             begin
               unten:=zneu[i].inf;
               oben:=zneu[i].sup;
               writeln('        [',unten:23:0:-1,',',oben:23:0:1,']');
             end;  
bis:= 0;
while (bis < N) and (zneu[bis+1].inf > 0) do bis:= bis + 1;
c:=bis*<schritt.inf;
if c > left then left:=c;
writeln('The free boundary c is greater than ',left:23:0:-1);
writeln('The free boundary c is less than ',a:23:0:1);
writeln(zaehler,' iteration steps have been done');
writeln('There were still ',count,' nondegenerate intervals');
                              
{ jetzt wollen wir mal sehen, ob wir a verbessern koennen }

bis:=N+1;
while (bis > 1) and (zneu[bis-1].sup =0) do bis:= bis - 1;
if bis < N+1 then a:=bis*>schritt.sup;
writeln('The free boundary c is less than ',a:23:0:1);           

end;


var N,k           : integer;
    y0,hilf,kappa : interval;
    a,sicherheit  : real;
    left          : real;

begin



hilf:=10;
y0:=1/hilf;
writeln('First we do the symmetric single step method');
writeln('with n=9');
writeln('Type in k');
writeln('then you consider N=99+k*100 components');
read(k);
N:=99+k*100;

kappa:=f_kappa(y0);
hilf:=sqrt(2*y0/kappa);
a:=hilf.sup;
left:=0;
repeat
sicherheit:=a;
mix(N,k,y0,left,a);
until sicherheit=a;
         
end.