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

Узнай цену своей работы

Формулировка задачи:

Доброго времени суток! Задание в прикрепленном файле. Есть уравнение: Необходимо вычислить
α
с помощью программ ZEROIN и QUANC8. Оформила подынтегральную функцию как function (starting_fun), затем написала еще одну функцию, в которой вызывается QUANC8 (fun). Ее значением должно быть значение интеграла минус единица. В основном теле программы вызываю ZEROIN (для промежутка [0;1]) со ссылкой на функцию с QUANC8 (т.е. fun), должно получаться то самое значение alpha, которое нужно найти, но.. QUANC8 не считает значения, выводит NAN. И все решение летит к чертям. При этом, когда в function подынтегральной функции (starting_fun) заменяю alpha на число, все в порядке (но, естественно, задача состоит именно в нахождении alpha). Получается, QUANC8 не видит этой переменной, я правильно понимаю? Как правильно сделать вызов этих функций? Или что можно еще сделать, чтобы считалось значение? В интернете пыталась искать что-то похожее, не получилось. Очень надеюсь на вашу помощь. Краткое описание ZEROIN
Спойлер
Краткое описание QUANC8
Спойлер
Заранее всем огромное спасибо!

Решение задачи: «Решение интегрального уравнения. QUANC8 и ZEROIN»

textual
Листинг программы
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.

Объяснение кода листинга программы

The code you provided is a Pascal program that implements a method for solving integral equations. It uses a user-defined function to evaluate the equation at different points. The program then adjusts the step size and repeats the process until the desired accuracy is achieved. The user-defined function UserFunction1 is used to evaluate the equation at a given point x. The user-defined function UserFunction3 is used to evaluate the equation at a different point x+d, where d is the step size. The program starts by initializing some variables and constants, including the step size step, the tolerance tolerr, and the number of function calls nofun. It then proceeds to the main loop, which iterates over different steps in the integration process. In each iteration, it calculates the current position x, evaluates the equation at x using UserFunction1, and updates the step size step. If the current step size step is smaller than the tolerance tolerr, the program increases the step size step by a factor of 2. Otherwise, it keeps the current step size. After each iteration, the program updates the variable result, which represents the sum of the current step size step and the tolerance tolerr. It also updates the variable errest, which represents the residual error between the calculated and exact values. The program continues iterating until the desired accuracy is achieved, which is determined by the tolerance tolerr. At the end of the program, the variable zeroin is set to the value of the final step size step. Please note that the code provided is quite complex and may require further explanation or adjustments to work correctly depending on the specific problem being solved.

Оцени полезность:

13   голосов , оценка 4.308 из 5
Похожие ответы