```(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *) (* Version 1, Update 2 *) MODULE CalcGrunwald; (** AUTHOR "adf"; PURPOSE "Grünwald-Letnikov algorithms for the fractional calculus"; *) IMPORT NbrInt, NbrRe, MathRe, CalcFn; CONST (** Status of integration, i.e., admissible values for the returned parameter 'res'. *) OKay* = 0; MaxSubDivReached* = 1; Oscillation* = 2; (** Arguments of differ-integration are: f the function being differintegrated, x the argument of f, also the upper limit of differintegration, order (> 0) the fractional order of differintegration, i.e., a, error the requested and achieved error tolerances, res status of the differintegration. *) VAR (** Upper bound for the number of subintervals for integration/differentiation. *) MaxIntervals-: NbrInt.Integer; PROCEDURE Grunwald( f: CalcFn.ReArg; a, b, order: NbrRe.Real; VAR error, result: NbrRe.Real; VAR res: NbrInt.Integer ); (* Use Grünwald-Letnikov definition for fractional-order differintegration. Requires a < b. f the argument of differintegration, a the lower limit of differintegration, b the upper limit of differintegration, order the fractional order of differintegration, i.e., a, if order > 0 then it is a fractional differentiation if order < 0 then it is a fractional integration error the requested and achieved error tolerances. *) VAR i, im1, ip1, n: NbrInt.Integer; absError, absX, dx, errorLast, fn, fnm1, last, tol: NbrRe.Real; fns, save: POINTER TO ARRAY OF NbrRe.Real; BEGIN (* Implements the G2 algorithm of K. B. Oldham and J. Spanier, "The Fractional Calculus", 1974, pg. 138. *) result := 0; n := 4; tol := NbrRe.Max( error, 100 * NbrRe.Epsilon ); NEW( fns, 3 ); dx := (b - a) / 2; FOR i := 0 TO 2 DO fns[i] := f( b - i * dx ) END; error := NbrRe.MaxNbr; LOOP save := fns; fns := NIL; NEW( fns, n + 1 ); dx := (b - a) / n; FOR i := 0 TO n DO IF NbrInt.Odd( i ) THEN fns[i] := f( b - i * dx ) ELSE fns[i] := save[i DIV 2] END END; last := result; result := 0; FOR i := n - 1 TO 1 BY -1 DO im1 := i - 1; ip1 := i + 1; (* G2 algorithm, pg. 138 in Oldham & Spanier. *) fn := order * (fns[ip1] - 2 * fns[i] + fns[im1]) / 2; fn := order * ((fns[ip1] - fns[im1]) - fn) / 4; fn := fns[i] - fn; result := (result + fn) * (im1 - order) / i END; fnm1 := f( b + dx ); fn := order * (fns[1] - 2 * fns[0] + fnm1) / 2; fn := order * ((fns[1] - fnm1) - fn) / 4; result := (result + fns[0] - fn) * MathRe.Power( n / b, order ); n := 2 * n; errorLast := error; absError := NbrRe.Abs( result - last ); absX := NbrRe.Abs( result ); IF absX > 1 THEN error := absError / absX ELSE error := NbrRe.Max( absError, 50 * NbrRe.Epsilon * absX ) END; IF error < tol THEN res := OKay; EXIT END; IF n = MaxIntervals THEN res := MaxSubDivReached; EXIT END; IF error > errorLast THEN result := last; res := Oscillation; EXIT END END END Grunwald; (** Computes a Riemann-Liouville fractional-order derivative. Daf(x) = DnIaf(x) = (1/G(n-a)) dn/dxn x0x (x-y)n-1-a f(y) dy, n-1 < a # n, n N @, a N !+ *) PROCEDURE SolveD*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR error: NbrRe.Real; VAR res: NbrInt.Integer ): NbrRe.Real; VAR result, zero: NbrRe.Real; BEGIN error := 0; result := 0; res := OKay; zero := 0; IF f # NIL THEN IF zero # x THEN IF order < zero THEN order := -order END; IF zero < x THEN Grunwald( f, zero, x, order, error, result, res ) ELSE Grunwald( f, x, zero, order, error, result, res ); result := -result END END ELSE result := 0 END; RETURN result END SolveD; (** Computes a Riemann-Liouville fractional-order integral. Iaf(x) = (1/G(a)) x0x (x-y)a-1 f(y) dy, a N !+ *) PROCEDURE SolveI*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR error: NbrRe.Real; VAR res: NbrInt.Integer ): NbrRe.Real; VAR result, zero: NbrRe.Real; BEGIN error := 0; result := 0; res := OKay; zero := 0; IF f # NIL THEN IF zero # x THEN IF order > zero THEN order := -order END; IF zero < x THEN Grunwald( f, zero, x, order, error, result, res ) ELSE Grunwald( f, x, zero, order, error, result, res ); result := -result END END END; RETURN result END SolveI; BEGIN MaxIntervals := NbrInt.MaxNbr END CalcGrunwald.```