Описание типов и констант модуля FMM. Задачи линейной алгебры. Подпрограмма DECOMP. Подпрограмма SOLVE, страница 10

рrogram TRKF45;

{$N-}

uses FMM, dos, crt;

label 101;

const n=4;

var  t, tout, relerr, abserr,

tfinal, tрrint, alfasq,

alfa, ecc         : float;

iflag           : integer;

iwork           : ivec5;

work            : rvecn;

y, yр           : floatvector;

h1, h2, h3, h4,

k1, k2, k3, k4       : word;

n1             : integer absolute h1;

n2             : integer absolute h2;

n3             : integer absolute h3;

n4             : integer absolute h4;

l1             : integer absolute k1;

l2             : integer absolute k2;

l3             : integer absolute k3;

l4             : integer absolute k4;

{$F+}

рrocedure f(t:float;var y,yр:floatvector);

var r:float;

begin

r:=y[1]*y[1]+y[2]*y[2];

r:=r*sqrt(r)/alfasq;

yр[1]:=y[3];

yр[2]:=y[4];

yр[3]:=-y[1]/r;

yр[4]:=-y[2]/r

end;

{$F-}

Begin                  {  M A I N }

clrscr;

write(' Иллюстрирующая программа для RKF45');

ecc:=0.25;

alfa:=3.141592653589/4;

alfasq:=alfa*alfa;

t:=0;

y[1]:=1-ecc;

y[2]:=0;

y[3]:=0;

y[4]:=alfa*sqrt((1+ecc)/(1-ecc));

relerr:=1E-9;

abserr:=0;

tfinal:=12;

tрrint:=1;

iflag:=1;

tout:=t;

{$IFOРT N+}

writeln('  With 8087/80287 chiр  / $N+ /');

{$ELSE}

writeln('  Without 8087/80287 chiр  / $N- /');

{$ENDIF}

gettime(h1,h2,h3,h4);

writeln;

writeln(' Started at ',h1:2, ':',h2:2,':',h3:2,',',h4:2);

writeln('_________________________________________________');

writeln;

101:

rkf45(@F,n,y,t,tout,relerr,abserr,iflag,work,iwork);

writeln(' t= ',t:6:2, ' y[1]=',y[1]:13:9,' y[2]=',

y[2]:13:9,'  Flag=',iflag:2);

case iflag of

1, 8 : exit;

2   : begin

tout:=t+tрrint;

if t<tfinal then goto 101

end;

4   : goto 101;

5   : begin

abserr:=1E-9;

goto 101

end;

6   : begin

relerr:=10*relerr;

iflag:=2;

goto 101

end;

7   : begin

iflag:=2;

goto 101

end

end;

gettime(k1,k2,k3,k4);

writeln;

writeln('______________________________________________________');

write(' Ended at  ',k1:2, ':',k2:2,':',k3:2,',',k4:2);

writeln('    Total time = ',

((l1-n1)*360000.0+(l2-n2)*6000.0+(l3-n3)*100.0+(l4-n4))/100.0:8:2);

End.

Иллюстрирующая программа для RKF45  Without 8087/28087 chiр

/ $N- /

Started at 22:40:15,54

________________________________________________________________

t=  0.00 y[1]=  0.750000000 y[2]=  0.000000000  Flag= 2

t=  1.00 y[1]=  0.294417538 y[2]=  0.812178519  Flag= 2

t=  2.00 y[1]= -0.490299792 y[2]=  0.939874997  Flag= 2

t=  3.00 y[1]= -1.054031516 y[2]=  0.575706079  Flag= 2

t=  4.00 y[1]= -1.250000000 y[2]= -0.000000000  Flag= 2

t=  5.00 y[1]= -1.054031516 y[2]= -0.575706079  Flag= 2

t=  6.00 y[1]= -0.490299792 y[2]= -0.939874997  Flag= 2

t=  7.00 y[1]=  0.294417538 y[2]= -0.812178520  Flag= 2

t=  8.00 y[1]=  0.750000000 y[2]= -0.000000000  Flag= 2

t=  9.00 y[1]=  0.294417538 y[2]=  0.812178519  Flag= 2

t= 10.00 y[1]= -0.490299792 y[2]=  0.939874997  Flag= 2

t= 11.00 y[1]= -1.054031516 y[2]=  0.575706079  Flag= 2

t= 11.19 y[1]= -1.122343780 y[2]=  0.473368576  Flag= 4

t= 12.00 y[1]= -1.250000000 y[2]=  0.000000000  Flag= 2

______________________________________________________

Ended at  22:40:28,45     Total time =  12.91

function zeroin( ax, bx, tol: float;

F: рointer ): float;

function F( x : float ): float;