rkqc.txt

用于开发税务票据管理的软件

procedure RKQC(var Y:array of real;var DYDX:array of real; N:integer;
      var X:real; HTRY, EPS:real; YSCAL:array of real;var HDID, HEND:real);
label 1;
const
    ONE = 1;  SAFETY = 0.9;  ERRCON = 0.0006;  FCOR = 0.066667;
    PGROW = -0.2;    PSHRNK = -0.25;
var
    YTEMP, YSAV, DYSAV:array[0..10] of real;
    XSAV,H,HH,ERRMAX:real;  I:integer;
begin
    XSAV:=X;
    For I:=1 To N do
    begin
        YSAV[I]:=Y[I];
        DYSAV[I]:=DYDX[I];
    end; 
    H:=HTRY;
1:  HH:=0.5 * H;
    RK4(YSAV, DYSAV, N, XSAV, HH, YTEMP);
    X:=XSAV + HH;
    DERIVS(X, YTEMP, DYDX);
    RK4(YTEMP, DYDX, N, X, HH, Y);
    X:=XSAV + H;
    If X = XSAV Then
    begin
      ShowMessage('Stepsize not significant in RKQC.');
      Exit;
    end;
    RK4(YSAV, DYSAV, N, XSAV, H, YTEMP);
    ERRMAX:=0;
    For I:=1 To N do
    begin
      YTEMP[I]:=Y[I] - YTEMP[I];
      If ERRMAX < Abs(YTEMP[I] / YSCAL[I]) Then
         ERRMAX:=Abs(YTEMP[I] / YSCAL[I]);
    end;
    ERRMAX:=ERRMAX / EPS;
    If ERRMAX > ONE Then
    begin
      H:=SAFETY * H * Exp(PSHRNK * Ln(ERRMAX));
      goto 1;
    end
    Else
    begin
      HDID:=H;
      If ERRMAX > ERRCON Then
        HEND:=SAFETY * H * Exp(PGROW * Ln(ERRMAX))
      Else
        HEND:=4 * H;
    end;
    For I:=1 To N do
      Y[I]:=Y[I] + YTEMP[I] * FCOR;
end;