program SymmetrischesEinzelschritt;

{ Beispiel 1 }
{ 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];
      y_alt,y_neu,ystart,y_halb  : ivector[1..n];
      sum1,sum2,sum3             : 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 begin
	 sum2[1]:=0;
	 if zaehl > 1 then sum1[1]:=sum3[1]
         else sum1[1]:=y_alt[2];
      end
      else if i=n then begin
	 sum2[n]:=y_halb[n-1];
	 sum1[n]:=0;
      end
      else begin
	 if zaehl > 1 then sum1[i]:=sum3[i]
         else sum1[i]:=y_alt[i+1];
         sum2[i]:=y_halb[i-1];
      end;
      y_halb[i]:=imax((sum1[i] + sum2[i] -q[i])/2)**y_alt[i];
   end;

   {jetzt y_neu berechnen und verwende sum2[1..n]}

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


      y_neu[i]:=imax((sum2[i] + sum3[i] -q[i])/2)**y_halb[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('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.