Решение интегрального уравнения. QUANC8 и ZEROIN - Turbo Pascal
Формулировка задачи:
- program cours;
- {$N+}
- uses FMM, crt;
- Var a, b, z, tol : double;
- i, j : integer;
- ai, bi, relerr, abserr,
- result, errest, flag : double;
- nofun : longint;
- function starting_fun(x : double) : double; { подынтегральная функция }
- var alpha : double;
- begin
- starting_fun:=1/(sqrt(2*(alpha + x*x*x/3 - x)));
- end;
- function fun(alpha: double) : double; { функция с вызовом quanc8 }
- begin
- ai:=0;
- bi:=1;
- relerr:=1.0e-12;
- abserr:=0.0;
- quanc8(@starting_fun,ai,bi,abserr,relerr,result,errest,nofun,flag);
- fun := result - 1;
- inc(j);
- end;
- BEGIN
- clrscr;
- a:=0;
- b:=1;
- tol:=1e-12;
- writeln(' Tol Root Nofun');
- writeln;
- j:=0;
- writeln(tol,' ',zeroin(a,b,tol,@fun),' ',j);
- END.
Решение задачи: «Решение интегрального уравнения. QUANC8 и ZEROIN»
- 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.
ИИ поможет Вам:
- решить любую задачу по программированию
- объяснить код
- расставить комментарии в коде
- и т.д