program Endversion; { f(x,s,t) = s + Arctan(xy') + pi } use i_ari, mvi_ari, mv_ari; 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; procedure G(x:rvector;ix,q:ivector;var init:boolean;var slopes:ivector); label 1; var argminy,argmaxz : rvector[lbound(q)..ubound(q)]; i,j : integer; nenner,zaehler,links : interval; wesen, hilf : interval; begin{function} i:=1; while i 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 slopes[i]:=-0.5 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; 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); slopes[i]:=-0.5+wesen*0.5; 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); slopes[i]:=-0.5 + wesen*0.5; 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); slopes[i]:=-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; 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); slopes[i]:=-wesen*0.5; 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); slopes[i]:= -wesen*0.5; 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); slopes[i]:=-wesen*0.5; end; end {3.2eigentlich} else begin {3.3} wesen:=intval(0,1); slopes[i]:=-wesen*0.5; end;{3.3} end;{3.2anfang} end;{3.Schritt} end;{2.Schritt} 1: i:=i+1; end;{while} init:= true; end;{function} function igatri(g,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)]:= 1; for k:= lbound(b) to ubound(b)-1 do begin ltri[k+1]:= g[k+1]/utri[k]; utri[k+1]:= 1 - ltri[k+1]*g[k]; 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]-g[i]*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 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 haupt(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 : ivector[1..n]; mitte : rvector[1..n]; newton, Hminvector : ivector[1..n]; rechteseite : ivector[1..n]; count, zaehler, zaehl : integer; initialisiert : boolean; begin kappa:=f_kappa(y0); Konst:=f_Konst(y0); Lipp:=f_Lipp(y0,a); zneu:=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; zaehler:=0; slopes:=0; initialisiert:=false; repeat zaehler:= zaehler+1; zalt:=zneu; for i:= 1 to n do mitte[i]:=zalt[i].inf; G(mitte,zalt,q,initialisiert,slopes); count:=0; for i:= 1 to n do if not (diam(slopes[i])=0) then count:=count+1; 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.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)* 0) do bis:= bis + 1; c:=bis* 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 : integer; y0,hilf,kappa : interval; a,sicherheit : real; left : real; begin hilf:=10; y0:=1/hilf; writeln('How many 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; haupt(n,y0,left,a); until sicherheit=a; end.