program Endversion;

{ f(x,s,t) = sqrt(1+s^2) }

use i_ari, mvi_ari, mv_ari;


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

function f_kappa(y0:interval):interval;
begin
f_kappa:=1;
end;

function f_Konst(y0:interval):interval;
var hilf:interval;
begin 
hilf:=sqrt(1+sqr(y0));
f_Konst:=hilf.sup;
end;

function f_Lipp(y0:interval):interval;
begin
f_Lipp:=0;
end;


function f(x,s,t:interval):interval;
begin
f:=sqrt(1+sqr(s));
end; 

function f_x(x,s,t:interval): interval;
begin
f_x:=0;
end;

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


function f_t(x,s,t:interval):interval;
begin
f_t:=0;
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 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 alg_C(n:integer;y0:interval;var left:real;var a:real);
var

      

      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         : integer;
    
begin

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

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

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;


zneu:=zstart;

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('slope matrix 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 5 = 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;
unten:=schritt.inf;
oben:=schritt.sup;
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 les than ',a:23:0:1);           

end;


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

begin



writeln('What is the initial value y0 ?');
read(y0);
writeln('How many ( > 2 !) components do you want to consider?'); 
read(n);

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