code 33080;
boolean procedure MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
value HMIN,HMAX,EPS,XEND,N,STIFF;
boolean FIRST,AVAILABLE,STIFF;
integer N;
real X,XEND,HMIN,HMAX,EPS;
array Y,YMAX,SAVE,JACOBIAN;
procedure DERIV,OUT;
begin own boolean ADAMS,WITH JACOBIAN;
own integer M,SAME,KOLD;
own real XOLD,HOLD,A0,TOLUP,TOL,TOLDWN,TOLCONV;
boolean EVALUATE,EVALUATED,DECOMPOSE,DECOMPOSED,CONV;
integer I,J,L,K,KNEW,FAILS;
real H, CH, CHNEW,ERROR,DFI,C;
array A[0:5],DELTA,LAST DELTA,DF[1:N],JAC[1:N, 1:N],AUX[1:3];
integer array P[1:N];
real procedure NORM2(AI); real AI;
begin real S,A; S:= 1.0"-100;
for I:= 1 step 1 until N do
begin A:= AI/YMAX[I]; S:= S + A * A end;
NORM2:= S
end NORM2;
procedure RESET;
begin if CH < HMIN/HOLD then CH:= HMIN/HOLD else
if CH > HMAX/HOLD then CH:= HMAX/HOLD;
X:= XOLD; H:= HOLD * CH; C:= 1;
for J:= 0 step M until K*M do
begin for I:= 1 step 1 until N do
Y[J+I]:= SAVE[J+I] * C;
C:= C * CH
end;
DECOMPOSED:= false
procedure METHOD;
begin I:= -39;
if ADAMS then
begin for C:= 1,1,144,4,0,.5,1,.5,576,144,1,5/12,1,
.75,1/6,1436,576,4,.375,1,11/12,1/3,1/24,
2844,1436,1,251/720,1,25/24,35/72,
5/48,1/120,0,2844,0.1
do begin I:= I+ 1; SAVE[I]:= C end
end else
begin for C:= 1,1,9,4,0,2/3,1,1/3,36,20.25,1,6/11,
1,6/11,1/11,84.028,53.778,0.25,.48,1,.7,.2,.02,
156.25, 108.51, .027778, 120/274, 1, 225/274,
85/274, 15/274, 1/274, 0, 187.69, .0047361
do begin I:= I + 1; SAVE[I]:= C end
end
end METHOD;
procedure ORDER;
begin C:= EPS * EPS; J:= (K-1) * (K + 8)/2 - 38;
for I:= 0 step 1 until K do A[I]:= SAVE[I+J];
TOLUP := C * SAVE[J + K + 1];
TOL := C * SAVE[J + K + 2];
TOLDWN := C * SAVE[J + K + 3];
TOLCONV:= EPS/(2 * N * (K + 2));
A0:= A[0]; DECOMPOSE:= true;
end ORDER;
procedure EVALUATE JACOBIAN;
begin EVALUATE:= false;
DECOMPOSE:= EVALUATED:= true;
if AVAILABLE then else
begin real D; array FIXY,FIXDY,DY[1:N];
for I:= 1 step 1 until N do
FIXY[I]:= Y[I];
DERIV(FIXDY);
for J:= 1 step 1 until N do
begin D:= if EPS > ABS(FIXY[J])
then EPS * EPS
else EPS * ABS(FIXY[J]);
Y[J]:= Y[J] + D; DERIV(DY);
for I:= 1 step 1 until N do
JACOBIAN[I,J]:= (DY[I]-FIXDY[I])/D;
Y[J]:= FIXY[J]
end
end
end EVALUATE JACOBIAN;
procedure DECOMPOSE JACOBIAN;
begin DECOMPOSE:= false;
DECOMPOSED:= true; C:= -A0 * H;
for J:= 1 step 1 until N do
begin for I:= 1 step 1 until N do
JAC[I,J]:= JACOBIAN[I,J] * C;
JAC[J,J]:= JAC[J,J] + 1
end;
AUX[2]:=1.0"-12;
DEC(JAC,N,AUX,P)
end DECOMPOSE JACOBIAN;
procedure CALCULATE STEP AND ORDER;
begin real A1,A2,A3;
A1:= if K <= 1 then 0 else
0.75 * (TOLDWN/NORM2(Y[K*M+I])) ** (0.5/K);
A2:= 0.80 * (TOL/ERROR) ** (0.5/(K + 1));
A3:= if K >= 5 or FAILS ^= 0
then 0 else
0.70 * (TOLUP/NORM2(DELTA[I] - LAST DELTA[I])) **
(0.5/(K+2));
if A1 > A2 and A1 > A3 then
begin KNEW:= K-1; CHNEW:= A1 end else
if A2 > A3 then
begin KNEW:= K ; CHNEW:= A2 end else
begin KNEW:= K+1; CHNEW:= A3 end
end CALCULATE STEP AND ORDER;
if FIRST then
begin FIRST:= false; M:= N;
for I:= -1,-2,-3 do SAVE[I]:= 0;
OUT(0,0);
ADAMS:= not STIFF; WITH JACOBIAN:= not ADAMS;
if WITH JACOBIAN then EVALUATE JACOBIAN;
METHOD;
NEW START: K:= 1; SAME:= 2; ORDER; DERIV(DF);
H:= if not WITH JACOBIAN then HMIN else
SQRT(2 * EPS/SQRT(NORM2 (MATVEC(1,N,I,JACOBIAN,DF))));
if H > HMAX then H:= HMAX else
if H < HMIN then H:= HMIN;
XOLD:= X; HOLD:= H; KOLD:= K; CH:= 1;
for I:= 1 step 1 until N do
begin SAVE[I]:= Y[I]; SAVE[M+I]:= Y[M+I]:= DF[I] * H
end;
OUT(0,0)
end else
begin WITH JACOBIAN:= not ADAMS; CH:= 1;
K:=KOLD; RESET; ORDER;
DECOMPOSE:= WITH JACOBIAN
end;
for L:= 0 while X < XEND do
begin if X + H <= XEND then X:= X + H else
begin H:= XEND-X; X:= XEND; CH:= H/HOLD; C:= 1;
for J:= M step M until K*M do
begin C:= C* CH;
for I:= J+1 step 1 until J+N do
Y[I]:= Y[I] * C
end;
SAME:= if SAME<3 then 3 else SAME+1;
end;
comment PREDICTION;
for L:= 1 step 1 until N do
begin for I:= L step M until (K-1)*M+L do
for J:= (K-1)*M+L step -M until I do
Y[J]:= Y[J] + Y[J+M];
DELTA[L]:= 0
end; EVALUATED:= false;
comment CORRECTION AND ESTIMATION LOCAL ERROR;
for L:= 1,2,3 do
begin DERIV(DF);
for I:= 1 step 1 until N do
DF[I]:= DF[I] * H - Y[M+I];
if WITH JACOBIAN then
begin if EVALUATE then EVALUATE JACOBIAN;
if DECOMPOSE then DECOMPOSE JACOBIAN;
SOL(JAC,N,P,DF)
end;
CONV:= true;
for I:= 1 step 1 until N do
begin DFI:= DF[I];
Y[ I]:= Y[ I] + A0 * DFI;
Y[M+I]:= Y[M+I] + DFI;
DELTA[I]:= DELTA[I] + DFI;
CONV:= CONV and ABS(DFI) < TOLCONV * YMAX[I]
end;
if CONV then
begin ERROR:= NORM2(DELTA[I]);
goto CONVERGENCE
end
comment ACCEPTANCE OR REJECTION;
if not CONV then
begin if not WITH JACOBIAN then
begin EVALUATE:= WITH JACOBIAN:= SAME >= K
or H<1.1 * HMIN;
if not WITH JACOBIAN then CH:= CH/4;
end else
if not DECOMPOSED then DECOMPOSE:= true else
if not EVALUATED then EVALUATE := true else
if H > 1.1 * HMIN then CH:= CH/4 else
if ADAMS then goto TRY CURTISS else
begin SAVE[-1]:= 1; goto RETURN end;
RESET
end else CONVERGENCE:
if ERROR > TOL then
begin FAILS:= FAILS + 1;
if H > 1.1 * HMIN then
begin if FAILS > 2 then
begin if ADAMS then
begin ADAMS:= false; METHOD end;
KOLD:= 0; RESET; goto NEW START
end else
begin CALCULATE STEP AND ORDER;
if KNEW ^= K then
begin K:= KNEW; ORDER end;
CH:= CH * CHNEW; RESET
end
end else
begin if ADAMS then TRY CURTISS:
begin ADAMS:= false; METHOD
end else
if K = 1 then
begin comment VIOLATE EPS CRITERION;
C:= EPS * SQRT(ERROR/TOL);
if C > SAVE[-3] then SAVE[-3]:= C;
SAVE[-2]:= SAVE[-2] + 1;
SAME:= 4; goto ERROR TEST OK
end;
K:= KOLD:= 1; RESET; ORDER; SAME:= 2
end
end else ERROR TEST OK:
FAILS:= 0;
for I:= 1 step 1 until N do
begin C:= DELTA[I];
for L:= 2 step 1 until K do
Y[L*M+I]:= Y[L*M+I] + A[L] * C;
if ABS(Y[I]) > YMAX[I] then
YMAX[I]:= ABS(Y[I])
end;
SAME:= SAME-1;
if SAME= 1 then
begin for I:= 1 step 1 until N do
LAST DELTA[I]:= DELTA[I]
end else
if SAME= 0 then
begin CALCULATE STEP AND ORDER;
if CHNEW > 1.1 then
begin DECOMPOSED:= false;
if K ^= KNEW then
begin if KNEW > K then
begin for I:= 1 step 1
until N do Y[KNEW*M+I]
:= DELTA[I] * A[K]/KNEW
end;
K:= KNEW; ORDER
end;
SAME:= K+1;
if CHNEW * H > HMAX
then CHNEW:= HMAX/H;
H:= H * CHNEW; C:= 1;
for J:= M step M until K*M do
begin C:= C * CHNEW;
for I:= J+1 step 1 until
J+N do Y[I]:= Y[I] * C
end
end
else SAME:= 10
end;
if X ^= XEND then
begin XOLD:= X; HOLD:= H; KOLD:= K; CH:= 1;
for I:= K * M + N step -1 until 1 do
SAVE[I]:= Y[I];
OUT(H,K)
end
end CORRECTION AND ESTIMATION LOCAL ERROR;
end STEP;
RETURN: SAVE[0]:= if ADAMS then 0 else 1;
MULTISTEP:= SAVE[-1]= 0 and SAVE[-2]=0
end MULTISTEP;
eop