(* ::Package:: *) (* :Name: NumericalMath`PolynomialFit` *) (* :Title: Curve Fitting over Polynomial Vector Spaces in a Numerically Stable Way *) (* :Author: Jerry B. Keiper *) (* :Summary: This package does polynomial least squares fitting of data and avoids the numerical instability that often results if the basis functions are chosen to be simple powers of a variable. Shifted Chebyshev polynomials of the first kind are used as a basis for the vector space of polynomials in a single variable. The shift is chosen based on the distribution of the abscissas of the data. The result is returned as a pure function, partly for convenience and partly to discourage the user from trying to rewrite it as a linear combination of simple powers. Such "simplification" would destroy much of the numerical stability gained by this method. *) (* :Context: NumericalMath`PolynomialFit` *) (* :Package Version: 2.0 *) (* :Copyright: Copyright 1990-2007, Wolfram Research, Inc. *) (* :History: Originally written by Jerry B. Keiper, March 1989. Revised by Jerry B. Keiper, December 1990. *) (* :Keywords: Fit, least squared error, polynomial approximation *) (* :Source: G. Dahlquist & A. Bjorck, Numerical Methods, Prentice-Hall, Englewood Cliffs, NJ, 1974. *) (* :Warnings: None. *) (* :Mathematica Version: 2.0 *) (* :Limitations: Although much more stable than Fit, PolynomialFit is still affected by floating-point round-off. It allows the user to evaluate the result at a symbolic argument (say x) and expand to get a result expressed in terms of basis {1, x, x^2, ..., x^n}. Such abuse is strongly discouraged due to the numerical stability that results. *) Message[General::obspkg, "NumericalMath`PolynomialFit`"] BeginPackage["NumericalMath`PolynomialFit`"] FittingPolynomial::usage = "FittingPolynomial[<>, n] represents a polynomial of degree n. It is returned \ by PolynomialFit." PolynomialFit::usage = "PolynomialFit[data, n] gives the least squares polynomial fit to data \ much as Fit[data,{1,x,x^2,...,x^n}, x] would do, except that the calculation \ is done in a way that is numerically stable and the result is given as a \ FittingPolynomial that is numerically stable. Changing the form of the result \ is NOT advised for reasons of numerical stability." Begin["NumericalMath`PolynomialFit`Private`"] issueObsoleteFunMessage[fun_, context_] := (Message[fun::obspkgfn, fun, context]; ) PolynomialFit[data_List, n_Integer] := (issueObsoleteFunMessage[PolynomialFit,"NumericalMath`PolynomialFit`"]; Module[{len = Length[data], answer}, FittingPolynomial[answer, n] /; (answer = ChebyFit[data, len, Min[n, len-1]]; answer =!= $Failed) ]); FittingPolynomial[f_Function, n_][x_] := (issueObsoleteFunMessage[FittingPolynomial,"NumericalMath`PolynomialFit`"]; f[x]); Format[FittingPolynomial[f_Function, n_]] := SequenceForm["FittingPolynomial[<>, ", Round[n], "]"]; Cheby[n_, x_?NumberQ] := ChebyshevT[n, x]; chebycoeff[Times[coeff_, cheby_Cheby]] := {coeff, cheby[[1]]}; chebycoeff[0.] = {0,0}; PolynomialFit::baddata = "`1` is not a list of ordered pairs of numbers." ChebyFit[data_, len_, n_] := Module[{i, x, t, a, b, c, chebylist, temp, prec}, If[ n < 0, Return[$Failed]]; If[VectorQ[data], If[ n == 0, Return[Function @@ {Apply[Plus,data]/len}]]; t = (2x-len-1)/(len-1), (* else *) a = Transpose[data]; If[n == 0, Return[Function @@ {Apply[Plus, a[[2]]]/len}]]; a = a[[1]]; If[!VectorQ[a], Message[PolynomialFit::baddata, data]; Return[$Failed]]; b = Apply[Plus,a]/len; If[!NumberQ[b], Message[PolynomialFit::baddata, data]; Return[$Failed]]; a = 2 Sqrt[(a . a)/len - b b]; t = (x-b)/a ]; prec = Precision[data]; If[prec === Infinity, prec = MachinePrecision]; t = N[t, prec]; chebylist = Fit[data, Table[Cheby[i, t], {i, 0, n}], x]; If[Head[chebylist] =!= Plus, Return[$Failed]]; (* form a simple list of coefficients including zeros *) chebylist = Map[chebycoeff, Apply[List, chebylist]]; temp = Table[0, {n+1}]; Do[temp[[chebylist[[i,2]]+1]] = chebylist[[i,1]], {i,Length[chebylist]}]; Clear[a,b,c]; (* These variables will be used to smuggle information inside of Function[ ] which is HoldAll. *) t = t /. x -> #; (* Construct the algorithm for evaluating the polynomial, but don't evaluate anything. *) Function[Module[{ta = 0, tb = 0, tc, tt = a, clist = b}, (* use Clenshaw's algorithm *) Do[tc = 2 tt tb - ta + clist[[ti]]; ta = tb; tb = tc, {ti, Round[c], Round[2], Round[-1]}]; tt tb - ta + First[clist]]] /. {a -> t, b -> temp, c -> n+1} ]; End[] (* NumericalMath`PolynomialFit`Private` *) EndPackage[] (* NumericalMath`PolynomialFit` *)