program singlestep;

{ f(x,s,t) = sqrt(1+sqr(y')) }

use i_ari, mvi_ari, mv_ari;



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

function f_Konst(y0:interval):interval;
var hilf:interval;
begin
f_Konst:=1;
end;

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


function f(x,s,t:interval):interval;
begin
f:=sqrt(1 + sqr(t));
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:=0;
end;


function f_t(x,s,t:interval):interval;
begin
f_t:=t/sqrt(1+sqr(t));
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 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;

procedure haupt(n:integer;y0:interval;var a:real);
var
      i,j,k,h,l,zaehl         : integer;
      q                       : ivector[1..n];
      hvector                 : ivector[1..n];
      y_alt,y_neu,ystart      : ivector[1..n];
      y                       : ivector[1..n];
      SB                      : ivector[1..n];
      unten,oben,kra,c        : real;
      hilf,schritt,h1,h2      : interval;
      Fausw                   : interval;
      Fstrich                 : real;
      arg1,arg2,arg3,part     : interval;
      fx,fs,ft                : interval;
      bis                     : integer;
      stop                    : integer;
      kappa,Konst,Lipp        : interval;




begin

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



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;


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



   for i:= 1 to n do ystart[i]:= intval(0,y0.sup);
   y_neu:=ystart;


zaehl:=0;

repeat
     zaehl:= zaehl+1;
     y_alt := y_neu;
     for i:= 1 to n do
     begin if i=1 then y_neu[i]:= imax((y_alt[i+1]-q[i])/2)**y_alt[i]
     else if i=n then  y_neu[i]:= imax((y_alt[i-1]-q[i])/2)**y_alt[i]
     else y_neu[i]:=imax((y_alt[i+1]+y_alt[i-1]-q[i])/2)**y_alt[i];
     end;



     if zaehl mod 10 = 0 then {neues [q]}
     begin

         bis:= 0;
         while (bis < n) and (y_neu[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:=y_neu[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(y_neu[i+1].inf,y_neu[i-1].sup);
     hilf:=Lipp*y_neu[i-1] + Konst*(a-(i-1)*schritt);
     arg3:=intval(-hilf.sup,0);
     end
     else
     begin
     arg1:=intval(0,2*>schritt.sup);
     arg2:=intval(y_neu[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;
     end;
     q[1]:=q[1]-y0;



     arg1:=intval((bis-1)*<schritt.inf,a);
     arg2:=intval(0,y_neu[bis-1].sup);
     hilf:=Lipp*y_neu[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;

     end
     else writeln('noch keine Verbesserung moeglich');


     end;

 until y_neu=y_alt;




writeln('We had ',zaehl,' iteration steps');


for i:= 1 to n do
begin
 unten:=y_neu[i].inf;
 oben:=y_neu[i].sup;
 writeln('      [',unten:23:0:-1,',',oben:23:0:1,']');
end;

c:=bis*<schritt.inf;


{ jetzt wollen wir mal sehen, ob wir a verbessern koennen }

bis:=n+1;
while (bis > 1) and (y_neu[bis-1].sup =0) do bis:= bis - 1;
if bis < n+1 then a:=bis*>schritt.sup;
writeln;
writeln('We have c in [',c:23:0:-1,',',a:23:0:1,']');





end;

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

begin



y0:=1/10;
writeln('How many ( > 2 ) components do you want to consider?');
read(n);
kappa:=f_kappa(y0);
hilf:=sqrt(2*y0/kappa);
a:=hilf.sup;

repeat
  sicherheit:=a;
  haupt(n,y0,a);
until sicherheit=a;

         
end.