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 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_neu[i-1]-q[i])/2)**y_alt[i] else y_neu[i]:=imax((y_alt[i+1]+y_neu[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.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)* 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.