DEFINITION MODULE ApproxLib;
(*------------------------------------------------------------------------*)
(* Verschiedene Routinen um Polynme oder Rationale Funktionen an eine *)
(* Funktion oder einen Datensatz anzupassen. *)
(* Module for fitting polynomals and splines to data sets *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: ApproxLib.def,v 1.4 2018/04/09 20:48:43 mriedl Exp mriedl $ *)
FROM LMathLib IMPORT Funktion1;
PROCEDURE Ausgleichsgerade(VAR A,B : LONGREAL;
VAR dA,dB : LONGREAL;
VAR Xq,Yq : LONGREAL;
VAR R : LONGREAL;
dim : CARDINAL;
T : LONGREAL;
VAR X,Y : ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Berechnet eine Ausgleichsgerade durch X,Y nach der Methode *)
(* der kleinsten Fehlerquadratsumme. *)
(* *)
(* A,B : f(x) = Ax + B *)
(* *)
(* R : Korellationsfaktor. *)
(* T : Student-Faktor *)
(* Xq,Yq : Mittelwerte von X,Y. *)
(*----------------------------------------------------------------*)
CONST FitStraightLine = Ausgleichsgerade; (* For english users ... *)
(*----------------------------------------------------------------*)
(* Calculates the coefficients of a least square fitted straigt *)
(* through the data provided in X,Y *)
(* *)
(* A,B : f(x) = Ax + B *)
(* *)
(* R : Corellationsfactor. *)
(* T : Student-Faktor *)
(* Xq,Yq : Mean values of X,Y *)
(*----------------------------------------------------------------*)
PROCEDURE TschebKoeff( U,O : LONGREAL;
VAR TKoeff : ARRAY OF LONGREAL;
Grad : CARDINAL;
Funkt : Funktion1);
(*----------------------------------------------------------------*)
(* Berechnet die Tschebytscheffkoeffizienten TKoeff einer T-Ent- *)
(* wicklung des Grades Grad f"ur die Funktion Funkt im Intervall *)
(* [U,O]. *)
(*----------------------------------------------------------------*)
PROCEDURE PolyKoeff(VAR PKoeff : ARRAY OF LONGREAL; (* Polynomkoeffizienten *)
VAR TKoeff : ARRAY OF LONGREAL; (* Tschebytscheffkoeff. *)
U,O : LONGREAL; (* Intervall d. Polynom *)
Grad : CARDINAL);
(*----------------------------------------------------------------*)
(* Transformiert die Koeffizienten einer Tschebytscheffentwick- *)
(* lung im Intervall [-1,1] in die entsprechenden Polynomkoef- *)
(* fizienten f"ur das Intervall [U,O] mit *)
(* f(x) = \sum_{i=0}^Grad PKoeff_i x^i *)
(*----------------------------------------------------------------*)
PROCEDURE RatKoeff( Funkt : Funktion1;
U,O : LONGREAL;
VAR Ak,Bk : ARRAY OF LONGREAL;
n,m : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet f"ur das Intervall [U,O] der Eingabefunktion *)
(* Funkt(x) die Koeffizienten der rationalen Funktion *)
(* R^{n,m}(x) = {\sum_{k=0}^n A_k x^k *)
(* \over *)
(* 1.0 + \sum_{k=1}^m B_k x^k} *)
(* die die Funkt(x) im quadratischen Mittel bestapproximiert. *)
(*----------------------------------------------------------------*)
PROCEDURE FitTscheb(dim : CARDINAL;
Grad : CARDINAL;
VAR X,Y : ARRAY OF LONGREAL;
VAR TKoeff : ARRAY OF LONGREAL;
VAR U,O : LONGREAL;
VAR StdAbw : LONGREAL;
VAR IFehl : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet die Entwicklungkoeffizienten TKoeff einer *)
(* Tschebytscheffentwicklung nach der Methode der minimalen *)
(* Fehlerquadratsumme als Anpassung an die Datenpaare X,Y *)
(* *)
(* dim : Anzahl der Datenpaare (X,Y) *)
(* Grad : Grad der T-Entwicklung *)
(* U,O : Intervall der T-Entwicklung (R"uckgabeparameter !) *)
(* IFehl : Siehe PROCEDURRE PolyFit. *)
(*----------------------------------------------------------------*)
PROCEDURE FitPoly( dim : CARDINAL; (* Anzahl Datenpaare X,Y *)
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
Grad : CARDINAL; (* Grad des Polynoms *)
VAR PKoeff : ARRAY OF LONGREAL;(* Koeff. des Polynoms *)
VAR StdAbw : LONGREAL;
VAR IFehl : CARDINAL);
(*----------------------------------------------------------------*)
(* Fittet ein Polynom des Grades Grad an die Datenpaare (X,Y) *)
(* nach der Methode der minimalen Fehlerquadratsumme. *)
(* *)
(* IFehl : *)
(* 1 : Anzahl der Datenpunkte kleiner zwei. *)
(* 2 : Anzahl der Datenpunkte zu gro3. *)
(* 3 : Grad des Polynoms kleiner eins. *)
(* 4 : Grad des Polynoms gr"o\3er als die Anzahl der *)
(* Datenpunkte. *)
(* 5 : Fehlergleichungssystem zu gro\3. *)
(* 6 : Keine L"osung gefunden. *)
(*----------------------------------------------------------------*)
PROCEDURE FitRat(VAR X,Y : ARRAY OF LONGREAL; (* Datenpunkte *)
dim : CARDINAL; (* Anzahl Datenpunkte *)
VAR Ak,Bk : ARRAY OF LONGREAL;
n,m : CARDINAL;
VAR StdAbw : LONGREAL); (* Standartabweichung *)
(*----------------------------------------------------------------*)
(* Berechent die Koeffizenten einer rationalen Funktion *)
(* R^{n,m}(x) = {\sum_{k=0}^n A_k x^k *)
(* \over *)
(* 1.0 + \sum_{k=1}^m B_k x^k} *)
(* die die Datenpunkt X,Y im Sinne der minimierten Fehler- *)
(* quadratsumme approximiert. *)
(* R_i = Y_i - Q^{m,n}(x) (Residuenvektor) *)
(*----------------------------------------------------------------*)
PROCEDURE SplineKoeff(VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet die kubischen Splinkoeffizeienten der in X,Y ueber- *)
(* gebenen Wertepaare in den Feldern B,C und D. Die Y[i] werden *)
(* als der konstanter Teil des Spline-Polynoms behandelt, d.h. *)
(* *)
(* $ f(x) = Y_i + B_i x + C_i x^2 + D x^3 $ *)
(* *)
(* Der Parameter N gibt die Anzahl der Datenpunkte in X,Y an. *)
(* *)
(* Calculates the coefficients of the cubic spline interpolation *)
(* through the N data points provided in the arrays X and Y. *)
(* The coefficients are returned in arrays B,C and D. The values *)
(* of Y form the constant part of the polinominals of degree 3 *)
(* *)
(* Ref.: Forsythe, George; Malcolm, Michael; Moler, Cleve; *)
(* "Computer Methods for Mathematical Computations", *)
(* Prentice-Hall, Upper Saddle River, US (1977) *)
(*----------------------------------------------------------------*)
PROCEDURE EvalSpline( x : LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Auswerten einer Spline-Interpolation am Punkt x [X_1<=x<=X_N] *)
(* Die Koeffizienten B,C und D koennen z.B. mit der Prozedur *)
(* SplineKoeff berechnet werden. Die Y[i] werden als der *)
(* konstanter Teil des Spline-Polynoms behandelt. *)
(* *)
(* $ f(x) = Y_i + B_i x + C_i x^2 + D x^3 $ *)
(* *)
(* Evaluate the cubic spline interpolation at point x where *)
(* *)
(* X[i] <= x <= X[i+1] , i = 0..N-1 *)
(* *)
(* The spline coefficients B,C and D should have been calculated *)
(* by procedure "SplineKoeff". *)
(* *)
(* Ref.: See SplineKoeff *)
(*----------------------------------------------------------------*)
PROCEDURE EvalSplineDerivative( x : LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR Y : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR D : ARRAY OF LONGREAL;
N : CARDINAL;
VAR y : LONGREAL;
VAR dy : LONGREAL;
VAR d2y : LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Auswerten einer Spline-Interpolation am Punkt x [X_1<=x<=X_N] *)
(* Es wird in y der Interpolationswert zurueckgegeben, in dy die *)
(* erste Ableitung am Punkt x und in d2y die zweite Ableitung. *)
(* iFehl = 1 wenn x ausserhalb des Wertebereichs liegen sollte. *)
(* *)
(* Calculates the spline-interpolation as in function EvalSpline *)
(* and returns the value in y. In dy the first and in d2y the *)
(* second derivative at point x are returend as well. *)
(* iFehl = 1 if x is outside of the interval [X_1 ... X_N] *)
(*----------------------------------------------------------------*)
PROCEDURE Spline(VAR Xa,Ya : ARRAY OF LONGREAL; (* Eingabevektoren *)
N : CARDINAL; (* Anzahl der Eingabewerte *)
X1,Xm : LONGREAL;
VAR Xn,Yn : ARRAY OF LONGREAL; (* Ausgabevektor *)
M : CARDINAL); (* Anzahl Splinepunkte *)
(*----------------------------------------------------------------*)
(* Berechnet aus den Werten in Xa,Ya kubisch splineinterpolierte *)
(* Werte in Xn,Yn im Bereich X1,XM. *)
(* *)
(* Xa,Ya : Vorgegebene Datenpunte / original data points *)
(* N : Anzahl der Datenpunkte / number of data points *)
(* X1,Xm : Bereich, in dem die berechneten Spline-Punkte *)
(* liegen sollen / range for the new values *)
(* Xn,Yn : Splineinterpolierte Wertepaare / new data points *)
(* M : Anzahl der erzeugten Datenpunkte in den Feldern *)
(* Xn,Yn / number of values to be calculated in Xn,Yn *)
(* *)
(* Calculates M cubic spline interpolated values in the range *)
(* X1 .. Xm. These interpolated values fit to the N original data *)
(* pairs provided in Xa,Ya. This routine is a driver routine for *)
(* SplineKoeff and EvalSpline. *)
(*----------------------------------------------------------------*)
CONST MaxR = 32; (* Max degree of nominator polynominal A[0..MaxL-1] *)
MaxL = 32; (* Max degree of denominator polynominal B[0..MaxL-1] *)
PROCEDURE RatTsch1( l,r : CARDINAL;
ap,ep : LONGREAL;
eps : LONGREAL;
Fct : Funktion1;
W : Funktion1;
VAR A,B : ARRAY OF LONGREAL;
VAR etabar : LONGREAL;
VAR iFehl : INTEGER);
(* ---------------------------------------------------------------*)
(* RatTsch1 determines the coefficients A[0:l] and B[0:r] of the *)
(* best fit rational function *)
(* *)
(* ra(x) := (A[0]+A[1]*x + ... + A[l]*x**l) / *)
(* (B[0]+B[1]*x + ....+ B[r]*x**r) *)
(* *)
(* for the function Fct and the weight function W, which *)
(* minimizes the maximum modulus of the weighted error function *)
(* *)
(* (Fct(x) - ra(x)) * W(x) *)
(* *)
(* on the interval [ap,ep]. Input parameters are the integer *)
(* degree l of the numerator of ra, the integer degree r of the *)
(* denominator of ra, the real parameters ap and ep specifying *)
(* the left and right endpoint of the approximation interval, *)
(* respectively, the real procedure Fct giving the function to be *)
(* approximated and the weight function W, also of type real *)
(* procedure. Fct must be continuous, W continuous and positive *)
(* on [ap,ep]. The real input parameter eps should be 10^(-t+1) *)
(* if t is the number of decimal digits in the mantissa of the *)
(* arithmetic used. Output parameters are the real arrays *)
(* *)
(* A[0:l], B[0:r] (with B[0]= 1), *)
(* *)
(* the real etabar with the property *)
(* *)
(* abs(etabar) = max(abs(Fct(x) - ra(x))*W(x) : ap <= x <= ep). *)
(* *)
(* The integer iFehl specifies, whether the approximation was *)
(* successful (iFehl=0) or not (iFehl=1,2,3,4,5): if iFehl=1, *)
(* then it was impossible to obtain a sufficiently accurate ra(x) *)
(* within 10 iterations. The values iFehl=2,3 or 4, indicate the *)
(* presence of a pole of the approximant ra(x) within [ap,ep]. *)
(* If iFehl=5 then the rational interpolation failed *)
(* *)
(* ==> l,r : Degree of rational approximation *)
(* See also definition of constans MaxR and MaxL *)
(* ==> ap,ep : Anfangspunkt,Endpunkt (start and end point) *)
(* ==> eps : see text *)
(* ==> Fct : Function to be approximated *)
(* ==> W : Weight function, if unsure us w(x) = 1 *)
(* <== A,B : Computes coefficients of the rational approx. *)
(* <== etabar : see text *)
(* <== iFehl : Error indicator *)
(* 0 : Approximation was successful *)
(* 1 : The number of iterations admitted by the *)
(* limit ito dit not suffice to complete the *)
(* Remes-Algorithm. Increasing the number *)
(* ito might help. *)
(* 2 : A rational function with a pole has *)
(* occured - deteted by hitting a point in *)
(* the search procedure where error function *)
(* > 10**6 or rational function > 10**20 *)
(* 3 : A pole of a rational function had been *)
(* detected - denominator changes sign *)
(* within the interval *)
(* 4 : Case 1 and 2 or 3 are combined *)
(* 5 : Rational interpolation fails *)
(* *)
(* HINT: If iFehl = -1 on ENTRY the nuber of iterations and the *)
(* return code of iFehl are printed to the console at the *)
(* end of the routine. *)
(* *)
(* This procedure is a Modula-2 translation of the Algol60 *)
(* procedure rattsch1 - a Algol60 version is available as well *)
(* *)
(* Ref: Werner, H; Stoer, J.; Bommas, W.; "Rational Chebyshev *)
(* Approximation", Handbook Series Methods of Approximation, *)
(* Numerische Mathematik, Springer, Berlin (DE) (1967) *)
(* *)
(* For L2 approximations the procedure ApproxLib.RatKoeff can *)
(* be used. *)
(* ---------------------------------------------------------------*)
END ApproxLib.