Решение интегрального уравнения. QUANC8 и ZEROIN - Turbo Pascal

UNIT FMM;
{$N+}
INTERFACE
 
 Const    ndim  = 21;{FOR LAB1=21; LAB2= ; LAB3= }
          ndimS = 10;
          ndimC = 10;
       maxnfe   = 3000; {„®ЇгбвЁ¬®Ґ зЁб«® ўлзЁб«Ґ*Ё© Yp=F(t,Y) ¤«п RKF45}
 
{$IFOPT N+}
type
  float = double; { 8 byte real, requires 8087 math chip}
{$ELSE}
type
  float = real;   { 6 byte real, no math chip required }
{$ENDIF}
 
type  floatvector   = array [1..ndim]           of float;
      floatmatrix   = array [1..ndim,1..ndim]   of float;
      rvecn         = array [1..6*ndim+3]       of float;
      ivec5         = array [1..5]              of integer;
      ivector       = array [1..ndim]           of integer;
      floatvectorS  = array [1..ndimS]          of float;
      floatvectorC  = array [1..ndimC]          of float;
      floatmatrixCS = array [1..ndimS,1..ndimC] of float;
 
 
 
  Procedure quanc8(  fun    : pointer;
          a,b,abserr,relerr : float;
              var result,
                  errest    : float;
              var nofun     : longint;
              var flag      : float);
 
 
function zeroin(       ax,bx,tol : float;
                       F         : pointer) : float;
 
function sign(     a,b    : float)   : float;
 
  IMPLEMENTATION
 
Type   Floatworkrecord   = record
        yp           : floatvector;
        h            : float;
        f1           : floatvector;
        f2           : floatvector;
        f3           : floatvector;
        f4           : floatvector;
        f5           : floatvector;
        savre        : float;
        savae        : float;
                   end;
      Integerworkrecord  = record
        nfe, kop, init, jflag, kflag   : integer;
                           end;
 
{$F+}
{$L Initval.OBJ}
function UserFunction1(     X : float;
                     ProcAddr : Pointer) : float; external;
 
Function UserFunction3(     X : float;
                     ProcAddr : Pointer) : float; external;
 
{$F-}
 
function sign(a,b:float):float;
begin
  if b=0 then sign:=0 else sign:=b/abs(b)*abs(a)
end;
 
Procedure quanc8(fun:pointer;a,b,abserr,relerr:float;var result,
                  errest:float;var nofun:longint;var flag:float);
label 72,90,30,50,62,60,70,75,80,82;
var
 w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp : float;
 qprev,qnow,qdiff,qleft,esterr,tolerr            : float;
 f,x                                             : array [1..16] of float;
 fsave,xsave                                     : array [1..8,1..30] of float;
 qright                                          : array [1..31] of float;
 levmin,levmax,levout,nomax,nofin,lev,nim,i,j,
 kk,stepp                                        : longint;
begin
{writeln('Q8 start a,b,abserr,relerr,result,errest, nofun',
#13,#10,a,b,abserr,relerr,result,errest, nofun);}
 levmin:=1;
 levmax:=30;
 levout:=6;
 nomax:=5000;
 stepp:=1; for kk:=1 to levout+1 do stepp:=stepp*2;
 nofin:=nomax-8*(levmax-levout+stepp);
 w0:=3956.0/14175.0;
 w1:=23552.0/14175.0;
 w2:=-3712.0/14175.0;
 w3:=41984.0/14175.0;
 w4:=-18160.0/14175.0;
 flag:=0.0;
 result:=0.0;
 cor11:=0.0;
 errest:=0.0;
 area:=0.0;
 nofun:=0;
 if a=b then goto 90;{end procedure}
 lev:=0;
 nim:=1;
 x0:=a;
 x[16]:=b;
 qprev:=0.0;
 f0:=userfunction1(x0,fun);
 stone:=(b-a)/16.0;
 x[8]:=(x0+x[16])/2.0;
 x[4]:=(x0+x[8])/2.0;
 x[12]:=(x[8]+x[16])/2.0;
 x[2]:=(x0+x[4])/2.0;
 x[6]:=(x[4]+x[8])/2.0;
 x[10]:=(x[8]+x[12])/2.0;
 x[14]:=(x[12]+x[16])/2.0;
 j:=2;
 while j<=16 do begin
                 f[j]:=userfunction1(x[j],fun);
                 j:=j+2;
                end;
 nofun:=9;
30:  x[1]:=(x0+x[2])/2.0;
     f[1]:=userfunction1(x[1],fun);
     j:=3;
     while j<=15 do begin
 
                     x[j]:=(x[j-1]+x[j+1])/2.0;
                     f[j]:=userfunction1(x[j],fun);
                     j:=j+2;
                    end;
 nofun:=nofun+8;
 step:=(x[16]-x0)/16.0;
 qleft:=(w0*(f0+f[8])+w1*(f[1]+f[7])+w2*(f[2]+f[6])+w3*(f[3]+f[5])+w4*f[4])*step;
 qright[lev+1]:=(w0*(f[8]+f[16])+w1*(f[9]+f[15])+w2*(f[10]+f[14])+w3*(f[11]+f[13])+w4*f[12])*step;
 qnow:=qleft+qright[lev+1];
 qdiff:=qnow-qprev;
 area:=area+qdiff;
 esterr:=abs(qdiff)/1023.0;
 if abserr>=(relerr*abs(area)) then tolerr:=abserr*(step/stone)
                             else tolerr:=relerr*abs(area)*(step/stone);
 if lev<levmin then goto 50;
 if lev>=levmax then goto 62;
 if nofun>nofin then goto 60;
 if esterr<=tolerr then goto 70;
50: nim:=2*nim;
    lev:=lev+1;
    for i:=1 to 8 do begin
                      fsave[i,lev]:=f[i+8];
                      xsave[i,lev]:=x[i+8];
                     end;
 qprev:=qleft;
 for i:=-1 downto -8 do begin
                   f[2*i+18]:=f[i+9];
                   x[2*i+18]:=x[i+9];
                  end;
 goto 30;
60: nofin:=2*nofin;
    levmax:=levout;
    flag:=flag+(b-x0)/(b-a);
    goto 70;
62: flag:=flag+1.0;
70: result:=result+qnow;
    errest:=errest+esterr;
    cor11:=cor11+qdiff/1023.0;
72: if nim=2*(nim div 2) then goto 75;
               nim:=nim div 2;
               lev:=lev-1;
    goto 72;
75: nim:=nim+1;
if lev<=0 then goto 80;
 qprev:=qright[lev];
 x0:=x[16];
 f0:=f[16];
 for i:=1 to 8 do begin
                   f[2*i]:=fsave[i,lev];
                   x[2*i]:=xsave[i,lev];
                  end;
goto 30;
80: result:=result+cor11;
    if errest=0.0 then goto 90;{end procedure}
82: temp:=abs(result)+errest;
    if temp<>abs(result) then goto 90;
    errest:=2.0*errest;
    goto 82;
90: exit; end;
 
function zeroin(ax,bx,tol : float; F : pointer) : float;
 
Label 10, 20, 30, 40, 50, 60, 70, 80, 90;
Const fierstcall : boolean = true;
      eps        : float   = 1.0;
var      a,b,c,d,e,fa,fb,fc,tol1,xm,p,q,r,s : float;
Begin
    if fierstcall then begin
10:                      eps:=eps/2.0;
                         tol1:=1.0+eps;
                         if(tol1 > 1.0)  then goto 10;
                       end;
    a:=ax;
    b:=bx;
    fa:=userfunction3(a,F);
    fb:=userfunction3(b,F);
20: c:=a;
    fc:=fa;
    d:=b-a;
    e:=d;
30: if(abs(fc) >= abs(fb)) then goto 40;
    a:=b;
    b:=c;
    c:=a;
    fa:=fb;
    fb:=fc;
    fc:=fa;
40: tol1:=2.0*eps*abs(b)+0.5*tol;
    xm:=0.5*(c-b);
    if(abs(xm) <= tol1)  then goto 90;
    if(fb = 0.0)  then goto 90;
    if(abs(e) < tol1)  then goto 70;
    if(abs(fa) <= abs(fb))  then goto 70;
    if(a <> c)  then goto 50;
    s:=fb/fa;
    p:=2.0*xm*s;
    q:=1.0-s;
    goto 60;
50: q:=fa/fc;
    r:=fb/fc;
    s:=fb/fa;
    p:=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
    q:=(q-1.0)*(r-1.0)*(s-1.0);
60: if(p > 0.0) then q:=-q;
    p:=abs(p);
    if((2.0*p) >= (3.0*xm*q-abs(tol1*q)))  then goto 70;
    if(p >= abs(0.5*e*q))  then goto 70;
    e:=d;
    d:=p/q;
    goto 80;
70: d:=xm;
    e:=d;
80: a:=b;
    fa:=fb;
    if(abs(d) > tol1)  then b:=b+d;
    if(abs(d) <= tol1)  then b:=b+sign(tol1,xm);
    fb:=userfunction3(b,F);
    if((fb*(fc/abs(fc))) > 0.0)  then goto 20;
    goto 30;
90: zeroin:=b;
end;
 
END.