(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *)
(* Version 1, Update 2 *)

MODULE CalcDiethelm;   (** AUTHOR "adf"; PURPOSE "Diethelm's algorithms for the fractional calculus"; *)

IMPORT NbrRe, DataErrors, MathRe, MathGamma, CalcFn, CalcD1;

(**  Arguments of differ-integration are:
		f(x)		the function being differ-integrated,
		x		the argument of f, also the upper limit of differ-integration,
		order	(> 0)  the fractional order of differ-integration, i.e., a,
		tol		the requested and achieved error tolerances.

	Arguments for fractional-order differential equations are:
		y		the dependent variable, i.e., the function whose solution is being sought,
		y0		the vector of initial conditions
					if  0  < a <= 1 then y0[0] =  y(0+)
					if  1 < a < 2 then y0[0] =  y(0+) and y0[1] =  y'(0+),
		f(x,y)	the forcing function, i.e., the right-hand side of the differential equation,
		x		the independent variable,
		order	(0 < order < 2)  the fractional order of differentiation, i.e., a,
		tol		the requested and achieved error tolerances.
	*)

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 );
		(* Insert a new row into the Romberg tableau and perform Richardson extrapolation on that row.
			If the tableau is full, then ratchet the first column up one row, throwing away the first entry,
			insert the newSolution into the last row, and perform Richardon extrapolation over the entire tablueau.  *)
		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
				(* The tableau is not full yet.  Add another row.  Build this row of the tableau. *)
				INC( index );  tableau[index, 0] := newSolution;
				IF index > 0 THEN RichardsonExtrapolation END
			ELSE
				(* The tableau is full.  Ratchet back the rows.  Insert the newSolution.  And rebuild the entire tableau. *)
				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;
			(* Update the error. *)
			IF index < 2 THEN error := NbrRe.MaxNbr
			ELSE
				(* Compute an absolute error. *)
				error := NbrRe.Abs( tableau[index, index - 1] - tableau[index - 1, index - 1] );
				IF NbrRe.Abs( tableau[index, index] ) > 1 THEN
					(* Compute its reletive error. *)
					error := NbrRe.Abs( error / tableau[index, index] )
				END
			END;
			(* Assign the enhanced solution. *)
			soln := tableau[index, index]
		END Update;

	END Romberg;

VAR
	minTol, maxTol: NbrRe.Real;

	PROCEDURE CreateDiffFactors( alpha: NbrRe.Real;  x: Romberg );
	(* Create the factors, i.e., 2^(r[k]), for Richardson extrapolation for differentiation.
		To be called after UpdateAWeights.  *)
	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
			(* These are guesses as to what they probably are.  Their actual values have not yet been derived. *)
			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 );
	(* Create the factors, i.e., 2^(r[k]), for Richardson extrapolation for integration.
		To be called after UpdateCWeights.  *)
	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 );
	(* Create the factors, i.e., 2^(r[k]), for Richardson extrapolation for differential equations.
		To be called after UpdateBWeights and UpdateCWeights.  *)
	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 (* tol is okay *) END
	END VerifyTolerance;

	(** Computes a Riemann-Liouville fractional-order integral.
		Iaf(x) = (1/G(a)) x0x (x-y)a-1 f(y) dy,  where  0 < a.  *)
	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  (* integrate *)
				CreateIntFactors( order, romberg );  gamma := MathGamma.Fn( 2 + order );
				(* The algorithmic loop. *)
				REPEAT
					UpdateCWeights;  n := LEN( cWeights^ ) - 1;  h := x / n;
					(* Function evaluations over the history. *)
					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;
					(* The quadrature algorithm. *)
					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;

	(** Computes a Caputo fractional-order derivative.
		D*af(x) = IaDnf(x) = (1/G(n-a)) x0x (x-y)n-1-a [dnf(y)/dyn] dy,  f(k)(0+) = f0[k],  where  0 < a < 2, a # 1. *)
	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  (* differentiate *)
				CreateDiffFactors( order, romberg );  gamma := MathGamma.Fn( 2 - order );
				(* Account for the initial conditions. *)
				f0 := f( NbrRe.Epsilon * x );
				IF order > 1 THEN df0 := CalcD1.Solve( f, 0, CalcD1.Forward ) ELSE df0 := 0 END;
				(* The algorithmic loop. *)
				REPEAT
					UpdateAWeights;  n := LEN( aWeights^ ) - 1;  h := x / n;  solution := 0;
					(* Function evaluations over the history. *)
					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;
					(* The quadrature algorithm. *)
					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;

	(** Solves a fractional-order differential equation of the Caputo type.
		D*ay(x) = f(x,y(x)),  y(k)(0+) = y0[k],  where  0  < a < 2. *)
	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  (* solve *)
				CreateDiffEqnFactors( order, romberg );  p := 1 + order;
				gamma1 := MathGamma.Fn( p );  gamma2 := MathGamma.Fn( 1 + p );
				(* The algorithmic loop. *)
				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
						(* Integrate along the solution path to get y(i*h) for i = 1...n. *)
						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.