code 33191;
procedure GMS(X, XE, R, Y, H, HMIN, HMAX, DELTA, DERIVATIVE,
JACOBIAN, AETA, RETA, N, JEV, LU, NSJEV,
LINEAR, OUT);
value R;
real X, XE, H, HMIN, HMAX, DELTA, AETA, RETA;
integer R, N, JEV, NSJEV, LU;
boolean LINEAR;
array Y;
procedure DERIVATIVE, JACOBIAN, OUT;
begin
integer I, J, K, L, NSJEV1, COUNT, COUNT1, KCHANGE;
real A, A1, ALFA, E, S1, S2, Z1, X0, XL0, XL1,
XL2, ETA, H0, H1, Q, Q1, Q2, Q12, Q22, Q1Q2, DISCR;
boolean UPDATE, CHANGE, REEVAL, STRATEGY;
integer array RI, CI[1:R];
array AUX[1:9], BD1, BD2[1:3,1:3], Y1,
Y0[1:R], HJAC, H2JAC2, RQZ[1:R,1:R], YL, FL[1:3 * R];
procedure INITIALIZATION;
begin LU:= JEV:= N:= NSJEV1:= KCHANGE:= 0; X0:= X; DISCR:= 0;
K:=1; H1:= H0:= H; COUNT:= -2; AUX[2]:= "-14; AUX[4]:= 8;
DUPVEC(1, R, 0, YL, Y); REEVAL:= CHANGE:= true;
STRATEGY:= HMIN ^= HMAX and ^LINEAR; Q1:= -1; Q2:= -2;
COUNT1:= 0; XL0:= XL1:= XL2:= 0
end INITIALIZATION;
procedure COEFFICIENT;
begin XL2:= XL1; XL1:= XL0; XL0:= X0;
if CHANGE then
begin if N > 2 then
begin Q1:= (XL1 - XL0) / H1;
Q2:= (XL2 - XL0) / H1
end;
Q12:= Q1 * Q1; Q22:= Q2 * Q2; Q1Q2:= Q1 * Q2;
A:= -(3 * ALFA + 1) / 12;
BD1[1,3]:= 1 + (1 / 3 - (Q1 + Q2) * .5) / Q1Q2;
BD1[2,3]:= (1 / 3 - Q2 * .5) / (Q12 - Q1Q2);
BD1[3,3]:= (1 / 3 - Q1 * .5) / (Q22 - Q1Q2);
BD2[1,3]:= -ALFA * .5 + A * (1 - Q1 - Q2) / Q1Q2;
BD2[2,3]:= A * (1 - Q2) / (Q12 - Q1Q2);
BD2[3,3]:= A * (1 - Q1) / (Q22 - Q1Q2);
if STRATEGY or N <= 2 then
begin BD1[2,2]:= 1 / (2 * Q1);
BD1[1,2]:= 1 - BD1[2,2];
BD2[2,2]:= -(3 * ALFA + 1) / (12 * Q1);
BD2[1,2]:= -BD2[2,2] - ALFA * .5
end
end
end COEFFICIENT;
procedure OPERATOR CONSTRUCTION;
begin if REEVAL then
begin JACOBIAN(HJAC, Y);
JEV:= JEV + 1; NSJEV1:= 0;
if DELTA <= -"15 then ALFA:= 1 / 3 else
begin Z1:= H1 * DELTA;
A:= Z1 * Z1 + 12; A1:= 6 * Z1;
if ABS(Z1) < .1 then
ALFA:= (Z1 * Z1 / 140 - 1) * Z1 / 30 else
if Z1 < -33 then
ALFA:= (A + A1) / (3 * Z1 * (2 + Z1)) else
begin E:= EXP(Z1); ALFA:= ((A - A1) *
E - A - A1) / (((2 - Z1) * E - 2 - Z1) *
Z1 * 3)
end
end;
S1:= -(1 + ALFA) * .5; S2:= (ALFA * 3 + 1) / 12
A:= H1 / H0; A1:= A * A;
if REEVAL then A:= H1;
if A ^= 1 then
for J:= 1 step 1 until R do
COLCST(1, R, J, HJAC, A);
for I:= 1 step 1 until R do
begin for J:= 1 step 1 until R do
begin Q:= H2JAC2[I,J]:= if REEVAL then
MATMAT(1, R, I, J, HJAC, HJAC)
else H2JAC2[I,J] * A1;
RQZ[I,J]:= S2 * Q
end;
RQZ[I,I]:= RQZ[I,I] + 1;
ELMROW(1, R, I, I, RQZ, HJAC, S1)
end;
GSSELM(RQZ, R, AUX, RI, CI); LU:= LU + 1;
REEVAL:= UPDATE:= false
end OPERATOR CONSTRUCTION;
procedure DIFFERENCE SCHEME;
begin if COUNT ^= 1 then
begin DUPVEC(1, R, 0, FL, YL);
DERIVATIVE(FL); N:= N + 1; NSJEV1:= NSJEV1 + 1
end;
MULVEC(1, R, 0, Y0, YL, (1 - ALFA) / 2 - BD1[1,K]);
for L:= 2 step 1 until K do
ELMVEC(1, R, R * (L - 1), Y0, YL, -BD1[L,K]);
for L:= 1 step 1 until K do
ELMVEC(1, R, R * (L - 1), Y0, FL, H1 * BD2[L,K]);
for I:= 1 step 1 until R do
Y[I]:= MATVEC(1, R, I, HJAC, Y0);
MULVEC(1, R, 0, Y0, YL, (1 - 3 * ALFA) / 12 - BD2[1,K]);
for L:= 2 step 1 until K do
ELMVEC(1, R, R * (L - 1), Y0, YL, -BD2[L,K]);
for I:= 1 step 1 until R do
Y[I]:= Y[I] + MATVEC(1, R, I, H2JAC2, Y0);
DUPVEC(1, R, 0, Y0, YL);
for L:= 1 step 1 until K do
ELMVEC(1, R, R * (L - 1), Y0, FL, H1 * BD1[L,K]);
ELMVEC(1, R, 0, Y, Y0, 1); SOLELM(RQZ, R, RI, CI, Y)
end DIFFERENCE SCHEME;
procedure NEXT INTEGRATION STEP;
begin for L:= 2, 1 do
begin DUPVEC(L * R + 1, (L + 1) * R, -R, YL, YL);
DUPVEC(L * R + 1, (L + 1) * R, -R, FL, FL)
end;
DUPVEC(1, R, 0, YL, Y)
end NEXT INTEGRATION STEP;
procedure TEST ACCURACY;
begin K:= 2;
DUPVEC(1, R, 0, Y1, Y); DIFFERENCE SCHEME; K:= 3;
ETA:= AETA + RETA * SQRT(VECVEC(1, R, 0, Y1, Y1));
ELMVEC(1, R, 0, Y, Y1, -1);
DISCR:= SQRT(VECVEC(1, R, 0, Y, Y));
DUPVEC(1, R, 0, Y, Y1)
end TEST ACCURACY;
procedure STEPSIZE;
begin X0:= X; H0:= H1;
if N <= 2 and ^LINEAR then K:= K + 1;
if COUNT = 1 then
begin A:= ETA / (.75 * (ETA + DISCR)) + .33;
H1:= if A <= .9 or A >= 1.1 then A * H0
else H0; COUNT:= 0;
REEVAL:= A <= .9 and NSJEV1 ^= 1;
COUNT1:= if A >= 1 or REEVAL then 0 else
COUNT1 + 1; if COUNT1 = 10 then
begin COUNT1:= 0; REEVAL:= true;
H1:= A * H0
end
end else
begin H1:= H; REEVAL:= NSJEV = NSJEV1 and
^STRATEGY and ^LINEAR
end;
if STRATEGY then H1:= if H1 > HMAX
then HMAX else if H1 < HMIN then HMIN else H1;
X:= X + H1; if X >= XE then
begin H1:= XE - X0; X:= XE end;
if N <= 2 and ^LINEAR then REEVAL:= true;
if H1 ^= H0 then
begin UPDATE:= true; KCHANGE:= 3 end;
if REEVAL then UPDATE:= true;
CHANGE:= KCHANGE > 0 and ^LINEAR;
KCHANGE:= KCHANGE - 1
end STEPSIZE;
INITIALIZATION; OUT; X:= X + H1;
OPERATOR CONSTRUCTION;
BD1[1,1]:= 1; BD2[1,1]:= -ALFA * .5;
if ^LINEAR then COEFFICIENT;
NEXT STEP: DIFFERENCE SCHEME;
if STRATEGY then COUNT:= COUNT + 1;
if COUNT = 1 then TEST ACCURACY;
OUT; if X >= XE then goto END;
STEPSIZE; if UPDATE then OPERATOR CONSTRUCTION;
if ^LINEAR then COEFFICIENT;
NEXT INTEGRATION STEP; goto NEXT STEP;
END:
end GMS;
eop