MODULE CalcDiethelm;
IMPORT NbrRe, DataErrors, MathRe, MathGamma, CalcFn, CalcD1;
TYPE
Romberg = OBJECT
VAR index: LONGINT;
error, soln: NbrRe.Real;
factors: POINTER TO ARRAY OF NbrRe.Real;
tableau: ARRAY 8 OF POINTER TO ARRAY OF NbrRe.Real;
PROCEDURE & Initialize*;
VAR i: LONGINT;
BEGIN
index := -1; error := 0; soln := 0; NEW( factors, 8 );
FOR i := 0 TO 7 DO NEW( tableau[i], i + 1 ) END
END Initialize;
PROCEDURE Update( newSolution: NbrRe.Real );
VAR i: LONGINT;
PROCEDURE RichardsonExtrapolation;
VAR k: LONGINT; term: NbrRe.Real;
BEGIN
FOR k := 1 TO index DO
term := tableau[index - 1, k - 1]; term := term - factors[k - 1] * tableau[index, k - 1];
term := term / (1 - factors[k - 1]); tableau[index, k] := term
END
END RichardsonExtrapolation;
BEGIN
IF index < 7 THEN
INC( index ); tableau[index, 0] := newSolution;
IF index > 0 THEN RichardsonExtrapolation END
ELSE
FOR i := 0 TO 6 DO tableau[i, 0] := tableau[i + 1, 0] END;
tableau[7, 0] := newSolution; index := 0;
FOR i := 1 TO 7 DO INC( index ); RichardsonExtrapolation END
END;
IF index < 2 THEN error := NbrRe.MaxNbr
ELSE
error := NbrRe.Abs( tableau[index, index - 1] - tableau[index - 1, index - 1] );
IF NbrRe.Abs( tableau[index, index] ) > 1 THEN
error := NbrRe.Abs( error / tableau[index, index] )
END
END;
soln := tableau[index, index]
END Update;
END Romberg;
VAR
minTol, maxTol: NbrRe.Real;
PROCEDURE CreateDiffFactors( alpha: NbrRe.Real; x: Romberg );
BEGIN
IF alpha < 1 THEN
x.factors[0] := MathRe.Power( 2, 2 - alpha ); x.factors[1] := 4;
x.factors[2] := MathRe.Power( 2, 3 - alpha ); x.factors[3] := MathRe.Power( 2, 4 - alpha );
x.factors[4] := 16; x.factors[5] := MathRe.Power( 2, 5 - alpha );
x.factors[6] := MathRe.Power( 2, 6 - alpha ); x.factors[7] := 64
ELSE
x.factors[0] := MathRe.Power( 2, 2 - alpha ); x.factors[1] := MathRe.Power( 2, 3 - alpha );
x.factors[2] := 4; x.factors[3] := MathRe.Power( 2, 4 - alpha );
x.factors[4] := MathRe.Power( 2, 5 - alpha ); x.factors[5] := 16;
x.factors[6] := MathRe.Power( 2, 6 - alpha ); x.factors[7] := MathRe.Power( 2, 7 - alpha )
END
END CreateDiffFactors;
PROCEDURE CreateIntFactors( alpha: NbrRe.Real; x: Romberg );
BEGIN
IF alpha < 1 THEN
x.factors[0] := 4; x.factors[1] := MathRe.Power( 2, 2 + alpha ); x.factors[2] := 8;
x.factors[3] := MathRe.Power( 2, 3 + alpha ); x.factors[4] := 16;
x.factors[5] := MathRe.Power( 2, 4 + alpha ); x.factors[6] := 32;
x.factors[7] := MathRe.Power( 2, 5 + alpha )
ELSIF alpha < 2 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := MathRe.Power( 2, 2 + alpha ); x.factors[3] := 16;
x.factors[4] := MathRe.Power( 2, 3 + alpha ); x.factors[5] := 32;
x.factors[6] := MathRe.Power( 2, 4 + alpha ); x.factors[7] := 64
ELSIF alpha < 3 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := MathRe.Power( 2, 2 + alpha );
x.factors[4] := 32; x.factors[5] := MathRe.Power( 2, 3 + alpha ); x.factors[6] := 64;
x.factors[7] := MathRe.Power( 2, 4 + alpha )
ELSIF alpha < 4 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := 32;
x.factors[4] := MathRe.Power( 2, 2 + alpha ); x.factors[5] := 64;
x.factors[6] := MathRe.Power( 2, 3 + alpha ); x.factors[7] := 128
ELSIF alpha < 5 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := 32; x.factors[4] := 64;
x.factors[5] := MathRe.Power( 2, 2 + alpha ); x.factors[6] := 128;
x.factors[7] := MathRe.Power( 2, 3 + alpha )
ELSIF alpha < 6 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := 32; x.factors[4] := 64;
x.factors[5] := 128; x.factors[6] := MathRe.Power( 2, 2 + alpha ); x.factors[7] := 256
ELSIF alpha < 7 THEN
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := 32; x.factors[4] := 64;
x.factors[5] := 128; x.factors[6] := 256; x.factors[7] := MathRe.Power( 2, 2 + alpha )
ELSE
x.factors[0] := 4; x.factors[1] := 8; x.factors[2] := 16; x.factors[3] := 32; x.factors[4] := 64;
x.factors[5] := 128; x.factors[6] := 256; x.factors[7] := 512
END
END CreateIntFactors;
PROCEDURE CreateDiffEqnFactors( alpha: NbrRe.Real; x: Romberg );
BEGIN
IF alpha < 1 THEN
x.factors[0] := MathRe.Power( 2, 1 + alpha ); x.factors[1] := 4;
x.factors[2] := MathRe.Power( 2, 2 + alpha ); x.factors[3] := MathRe.Power( 2, 3 + alpha );
x.factors[4] := 16; x.factors[5] := MathRe.Power( 2, 4 + alpha );
x.factors[6] := MathRe.Power( 2, 5 + alpha ); x.factors[7] := 64
ELSE
x.factors[0] := 4; x.factors[1] := MathRe.Power( 2, 1 + alpha );
x.factors[2] := MathRe.Power( 2, 2 + alpha ); x.factors[3] := 16;
x.factors[4] := MathRe.Power( 2, 3 + alpha ); x.factors[5] := MathRe.Power( 2, 4 + alpha );
x.factors[6] := 64; x.factors[7] := MathRe.Power( 2, 5 + alpha )
END
END CreateDiffEqnFactors;
PROCEDURE VerifyTolerance( VAR tol: NbrRe.Real );
BEGIN
tol := NbrRe.Abs( tol );
IF tol < minTol THEN tol := minTol ELSIF tol > maxTol THEN tol := maxTol ELSE END
END VerifyTolerance;
PROCEDURE SolveI*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR tol: NbrRe.Real ): NbrRe.Real;
VAR i, n: LONGINT; firstWeight, gamma, h, solution: NbrRe.Real;
cWeights, fn, save: POINTER TO ARRAY OF NbrRe.Real; pWeights: ARRAY 3 OF NbrRe.Real;
romberg: Romberg;
PROCEDURE UpdateCWeights;
VAR len, newLen: LONGINT; m: NbrRe.Real;
BEGIN
m := 1 + order;
IF cWeights = NIL THEN
NEW( cWeights, 5 );
pWeights[0] := MathRe.Power( 2, m );
pWeights[1] := MathRe.Power( 3, m );
pWeights[2] := MathRe.Power( 4, m );
cWeights[0] := 1;
cWeights[1] := pWeights[0] - 2;
cWeights[2] := pWeights[1] - 2 * pWeights[0] + 1;
cWeights[3] := pWeights[2] - 2 * pWeights[1] + pWeights[0];
pWeights[0] := pWeights[1];
pWeights[1] := pWeights[2];
pWeights[2] := MathRe.Power( 5, m );
cWeights[4] := pWeights[2] - 2 * pWeights[1] + pWeights[0];
firstWeight := (m / 4 - 1 ) * pWeights[1] + pWeights[0]
ELSE
save := cWeights; len := LEN( cWeights^ ); cWeights := NIL;
newLen := 2 * (len - 1) + 1; NEW( cWeights, newLen );
FOR i := 0 TO len - 1 DO cWeights[i] := save[i] END;
FOR i := len TO newLen - 1 DO
pWeights[0] := pWeights[1];
pWeights[1] := pWeights[2];
pWeights[2] := MathRe.Power( i + 1, m );
cWeights[i] := pWeights[2] - 2 * pWeights[1] + pWeights[0]
END;
firstWeight := (m / (newLen - 1) - 1) * pWeights[1] + pWeights[0]
END
END UpdateCWeights;
BEGIN
IF order > 0 THEN
NEW( romberg ); VerifyTolerance( tol );
IF x < 0 THEN solution := -SolveI( f, -x, order, tol )
ELSIF x = 0 THEN solution := 0
ELSE
CreateIntFactors( order, romberg ); gamma := MathGamma.Fn( 2 + order );
REPEAT
UpdateCWeights; n := LEN( cWeights^ ) - 1; h := x / n;
IF n = 4 THEN
NEW( fn, 5 );
FOR i := 0 TO n DO fn[i] := f( h * i ) END
ELSE
save := fn; fn := NIL; NEW( fn, n + 1 );
FOR i := 0 TO n BY 2 DO fn[i] := save[i DIV 2] END;
FOR i := 1 TO n-1 BY 2 DO fn[i] := f( h * i ) END;
END;
solution := firstWeight * fn[0];
FOR i := 1 TO n DO solution := solution + cWeights[n - i] * fn[i] END;
romberg.Update( MathRe.Power( h, order ) * solution / gamma )
UNTIL romberg.error < tol;
solution := romberg.soln; tol := romberg.error
END
ELSE solution := 0; DataErrors.ReError( order, "The requested order of integration is not allowed." )
END;
RETURN solution
END SolveI;
PROCEDURE SolveD*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR tol: NbrRe.Real ): NbrRe.Real;
VAR i, n: LONGINT; df0, f0, gamma, h, lastWeight, solution: NbrRe.Real;
aWeights, fn, save: POINTER TO ARRAY OF NbrRe.Real; mWeights: ARRAY 3 OF NbrRe.Real;
romberg: Romberg;
PROCEDURE UpdateAWeights;
VAR len, newLen: LONGINT; m: NbrRe.Real;
BEGIN
m := 1 - order;
IF aWeights = NIL THEN
NEW( aWeights, 5 );
mWeights[0] := MathRe.Power( 2, m );
mWeights[1] := MathRe.Power( 3, m );
mWeights[2] := MathRe.Power( 4, m );
aWeights[0] := 1;
aWeights[1] := mWeights[0] - 2;
aWeights[2] := mWeights[1] - 2 * mWeights[0] + 1;
aWeights[3] := mWeights[2] - 2 * mWeights[1] + mWeights[0];
mWeights[0] := mWeights[1];
mWeights[1] := mWeights[2];
mWeights[2] := MathRe.Power( 5, m );
aWeights[4] := mWeights[2] - 2 * mWeights[1] + mWeights[0];
lastWeight := (m / 4 - 1) * mWeights[1] + mWeights[0]
ELSE
save := aWeights; len := LEN( aWeights^ ); aWeights := NIL;
newLen := 2 * (len - 1) + 1; NEW( aWeights, newLen );
FOR i := 0 TO len - 1 DO aWeights[i] := save[i] END;
FOR i := len TO newLen - 1 DO
mWeights[0] := mWeights[1];
mWeights[1] := mWeights[2];
mWeights[2] := MathRe.Power( i + 1, m );
aWeights[i] := mWeights[2] - 2 * mWeights[1] + mWeights[0]
END;
lastWeight := (m / (newLen - 1) - 1) * mWeights[1] + mWeights[0]
END
END UpdateAWeights;
BEGIN
IF ((0 < order) & (order # 1) & (order < 2)) THEN
NEW( romberg ); VerifyTolerance( tol );
IF x < 0 THEN solution := -SolveD( f, -x, order, tol )
ELSIF x = 0 THEN solution := 0
ELSE
CreateDiffFactors( order, romberg ); gamma := MathGamma.Fn( 2 - order );
f0 := f( NbrRe.Epsilon * x );
IF order > 1 THEN df0 := CalcD1.Solve( f, 0, CalcD1.Forward ) ELSE df0 := 0 END;
REPEAT
UpdateAWeights; n := LEN( aWeights^ ) - 1; h := x / n; solution := 0;
IF n = 4 THEN
NEW( fn, 5 );
FOR i := 0 TO n DO fn[i] := f( h * (n - i) ) END
ELSE
save := fn; fn := NIL; NEW( fn, n + 1 );
FOR i := 0 TO n BY 2 DO fn[i] := save[i DIV 2] END;
FOR i := 1 TO n-1 BY 2 DO fn[i] := f( h * (n - i) ) END;
END;
FOR i := 0 TO n - 1 DO
solution := solution + aWeights[i] * (fn[i] - f0 - h * (n - i) * df0)
END;
solution := solution + lastWeight * (fn[n] - f0);
romberg.Update( solution / (gamma * MathRe.Power( h, order )) )
UNTIL romberg.error < tol;
solution := romberg.soln; tol := romberg.error
END
ELSE solution := 0; DataErrors.ReError( order, "The requested order of differentiation is not allowed." )
END;
RETURN solution
END SolveD;
PROCEDURE SolveFODE*( f: CalcFn.Re2Arg; y0: ARRAY OF NbrRe.Real; x, order: NbrRe.Real;
VAR tol: NbrRe.Real ): NbrRe.Real;
VAR i, k, n: LONGINT; cCoef, pCoef, cSum, pSum, firstWeight, fn, gamma1, gamma2, h,
h2alpha, initialCondition, p, predictor, solution: NbrRe.Real;
bWeights, cWeights, save, y: POINTER TO ARRAY OF NbrRe.Real;
weights: ARRAY 2 OF NbrRe.Real; pWeights: ARRAY 3 OF NbrRe.Real;
romberg: Romberg;
PROCEDURE UpdateBWeights;
VAR len, newLen: LONGINT; m: NbrRe.Real;
BEGIN
m := order;
IF bWeights = NIL THEN
NEW( bWeights, 5 );
weights[0] := MathRe.Power( 2, m );
weights[1] := MathRe.Power( 3, m );
bWeights[0] := 0;
bWeights[1] := 1;
bWeights[2] := weights[0] - 1;
bWeights[3] := weights[1] - weights[0];
weights[0] := weights[1];
weights[1] := MathRe.Power( 4, m );
bWeights[4] := weights[1] - weights[0]
ELSE
save := bWeights; len := LEN( bWeights^ ); bWeights := NIL;
newLen := 2 * (len - 1) + 1; NEW( bWeights, newLen );
FOR i := 0 TO len - 1 DO bWeights[i] := save[i] END;
FOR i := len TO newLen - 1 DO
weights[0] := weights[1];
weights[1] := MathRe.Power( i, m );
bWeights[i] := weights[1] - weights[0]
END
END
END UpdateBWeights;
PROCEDURE UpdateCWeights;
VAR len, newLen: LONGINT; m: NbrRe.Real;
BEGIN
m := 1 + order;
IF cWeights = NIL THEN
NEW( cWeights, 5 );
pWeights[0] := MathRe.Power( 2, m );
pWeights[1] := MathRe.Power( 3, m );
pWeights[2] := MathRe.Power( 4, m );
cWeights[0] := 1;
cWeights[1] := pWeights[0] - 2;
cWeights[2] := pWeights[1] - 2 * pWeights[0] + 1;
cWeights[3] := pWeights[2] - 2 * pWeights[1] + pWeights[0];
pWeights[0] := pWeights[1];
pWeights[1] := pWeights[2];
pWeights[2] := MathRe.Power( 5, m );
cWeights[4] := pWeights[2] - 2 * pWeights[1] + pWeights[0]
ELSE
save := cWeights; len := LEN( cWeights^ ); cWeights := NIL;
newLen := 2 * (len - 1) + 1; NEW( cWeights, newLen );
FOR i := 0 TO len - 1 DO cWeights[i] := save[i] END;
FOR i := len TO newLen - 1 DO
pWeights[0] := pWeights[1];
pWeights[1] := pWeights[2];
pWeights[2] := MathRe.Power( i + 1, m );
cWeights[i] := pWeights[2] - 2 * pWeights[1] + pWeights[0]
END
END
END UpdateCWeights;
BEGIN
IF ((0 < order) & (order < 2)) THEN
NEW( romberg ); VerifyTolerance( tol );
IF order <= 1 THEN
IF LEN( y0 ) < 1 THEN
DataErrors.Error( "Initial-condition vector is of insufficient length." ); RETURN 0
END
ELSE
IF LEN( y0 ) < 2 THEN
DataErrors.Error( "Initial-condition vector is of insufficient length." ); RETURN 0
END
END;
IF x < 0 THEN solution := -SolveFODE( f, y0, -x, order, tol )
ELSIF x = 0 THEN solution := y0[0]
ELSE
CreateDiffEqnFactors( order, romberg ); p := 1 + order;
gamma1 := MathGamma.Fn( p ); gamma2 := MathGamma.Fn( 1 + p );
REPEAT
UpdateBWeights; UpdateCWeights; n := LEN( cWeights^ ) - 1;
h := x / n; h2alpha := MathRe.Power( h, order );
pCoef := h2alpha / gamma1; cCoef := h2alpha / gamma2;
y := NIL; NEW( y, n + 1 ); y[0] := y0[0];
FOR i := 1 TO n DO
initialCondition := y0[0]; firstWeight := (p / i - 1) * MathRe.Power( i, p );
IF order > 1 THEN initialCondition := initialCondition + i * h * y0[1] END;
IF i > 1 THEN firstWeight := firstWeight + MathRe.Power( i - 1, p ) END;
fn := f( 0, y[0] ); pSum := bWeights[i] * fn; cSum := firstWeight * fn;
FOR k := 1 TO i - 1 DO
fn := f( k * h, y[k] );
pSum := pSum + bWeights[i - k] * fn; cSum := cSum + cWeights[i - k] * fn
END;
predictor := initialCondition + pCoef * pSum;
y[i] := initialCondition + cCoef * (cSum + cWeights[0]* f( i * h, predictor ))
END;
romberg.Update( y[n] )
UNTIL romberg.error < tol;
solution := romberg.soln; tol := romberg.error
END
ELSE solution := 0; DataErrors.ReError( order, "The requested order of differentiation is not allowed." )
END;
RETURN solution
END SolveFODE;
BEGIN
minTol := MathRe.Sqrt( NbrRe.Epsilon ); maxTol := 0.1
END CalcDiethelm.