MODULE CalcD2;
IMPORT NbrInt, NbrRe, NbrCplx, MathRe, CalcFn;
CONST
Forward* = 9; Central* = 10; Backward* = 11;
VAR
epsilon, zero: NbrRe.Real;
PROCEDURE DoNothing( x: NbrRe.Real );
END DoNothing;
PROCEDURE DoCplxNothing( z: NbrCplx.Complex );
END DoCplxNothing;
PROCEDURE Solve*( f: CalcFn.ReArg; atX: NbrRe.Real; differencing: NbrInt.Integer ): NbrRe.Real;
VAR h, h2, hOpt, hMin, power, result, temp: NbrRe.Real;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrRe.Abs( atX ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
temp := atX + h; DoNothing( temp ); h := temp - atX; h2 := h * h;
IF differencing = Forward THEN
result := f( atX + 2 * h );
result := result - 2 * f( atX + h );
result := (result + f( atX )) / h2
ELSIF differencing = Backward THEN
result := f( atX );
result := result - 2 * f( atX - h );
result := (result + f( atX - 2 * h )) / h2
ELSE
result := f( atX + h );
result := result - 2 * f( atX );
result := (result + f( atX - h )) / h2
END;
RETURN result
END Solve;
PROCEDURE SolveCplx*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex;
VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, result, temp: NbrCplx.Complex;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
NbrCplx.Set( h, h, ch );
temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch;
IF differencing = Forward THEN
result := f( atZ + 2 * ch );
result := result - 2 * f( atZ + ch );
result := (result + f( atZ )) / ch2
ELSIF differencing = Backward THEN
result := f( atZ );
result := result - 2 * f( atZ - ch );
result := (result + f( atZ - 2 * ch )) / ch2
ELSE
result := f( atZ + ch );
result := result - 2 * f( atZ );
result := (result + f( atZ - ch )) / ch2
END;
RETURN result
END SolveCplx;
PROCEDURE SolveCplxRe*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex;
VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, result, temp: NbrCplx.Complex;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
NbrCplx.Set( h, zero, ch );
temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch;
IF differencing = Forward THEN
result := f( atZ + 2 * ch );
result := result - 2 * f( atZ + ch );
result := (result + f( atZ )) / ch2
ELSIF differencing = Backward THEN
result := f( atZ );
result := result - 2 * f( atZ - ch );
result := (result + f( atZ - 2 * ch )) / ch2
ELSE
result := f( atZ + ch );
result := result - 2 * f( atZ );
result := (result + f( atZ - ch )) / ch2
END;
RETURN result
END SolveCplxRe;
PROCEDURE SolveCplxIm*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex;
VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, result, temp: NbrCplx.Complex;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
NbrCplx.Set( zero, h, ch );
temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch;
IF differencing = Forward THEN
result := f( atZ + 2 * ch );
result := result - 2 * f( atZ + ch );
result := (result + f( atZ )) / ch2
ELSIF differencing = Backward THEN
result := f( atZ );
result := result - 2 * f( atZ - ch );
result := (result + f( atZ - 2 * ch )) / ch2
ELSE
result := f( atZ + ch );
result := result - 2 * f( atZ );
result := (result + f( atZ - ch )) / ch2
END;
RETURN result
END SolveCplxIm;
PROCEDURE SolveCplxAbs*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex;
VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, result, temp: NbrCplx.Complex;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrCplx.Abs( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
NbrCplx.SetPolar( h, zero, ch );
temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch;
IF differencing = Forward THEN
result := f( atZ + 2 * ch );
result := result - 2 * f( atZ + ch );
result := (result + f( atZ )) / ch2
ELSIF differencing = Backward THEN
result := f( atZ );
result := result - 2 * f( atZ - ch );
result := (result + f( atZ - 2 * ch )) / ch2
ELSE
result := f( atZ + ch );
result := result - 2 * f( atZ );
result := (result + f( atZ - ch )) / ch2
END;
RETURN result
END SolveCplxAbs;
PROCEDURE SolveCplxArg*( f: CalcFn.CplxArg; atZ: NbrCplx.Complex; differencing: NbrInt.Integer ): NbrCplx.Complex;
VAR h, hOpt, hMin, power: NbrRe.Real; ch, ch2, result, temp: NbrCplx.Complex;
BEGIN
power := 4 / 5; hMin := MathRe.Power( NbrRe.Epsilon, power ); power := 1 / 4;
hOpt := NbrCplx.Arg( atZ ) * MathRe.Power( epsilon, power ); h := NbrRe.Max( hOpt, hMin );
NbrCplx.SetPolar( zero, h, ch );
temp := atZ + ch; DoCplxNothing( temp ); ch := temp - atZ; ch2 := ch * ch;
IF differencing = Forward THEN
result := f( atZ + 2 * ch );
result := result - 2 * f( atZ + ch );
result := (result + f( atZ )) / ch2
ELSIF differencing = Backward THEN
result := f( atZ );
result := result - 2 * f( atZ - ch );
result := (result + f( atZ - 2 * ch )) / ch2
ELSE
result := f( atZ + ch );
result := result - 2 * f( atZ );
result := (result + f( atZ - ch )) / ch2
END;
RETURN result
END SolveCplxArg;
BEGIN
epsilon := 100 * NbrRe.Epsilon; zero := 0
END CalcD2.