12.   Curve Fitting  and Interpolation

2015-05-22

$Version

Often it is necessary to find mathematical expressions for data given as discrete points  only:

				{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])}, ..., {(Subscript[x, s]),(Subscript[y, s])}}.
				
 There are several methods available for performing this task. One may either choose to have 
 1.	one expression defined  over the whole interval containing all the points: This is done by fits or interpolating  polynomials;     	 
 	or  
 2.	one may proceed  such that different expressions are sought for the different intervals between the points;  the objects so 
 	obtained are called   interpolating function objects  or  splines.

12.1  Fitting Functions to Given Data

A function f(x)  is defined as a linear sum of given functions   (Subscript[ϕ, i])(x)   with unknown coefficients  (Subscript[c, i]) :

			f(x)  :=   (Subscript[c, 1]) (Subscript[ϕ, 1])(x)  +  (Subscript[c, 2]) (Subscript[ϕ, 2])(x)  +  ...  +  (Subscript[c, n]) (Subscript[ϕ, n])(x) .   			(12.1)
			
The number of functions, n,  must not  be larger than the number of data pairs, s. In many cases it is be smaller. Setting (Subscript[y, i]) =  f((Subscript[x, i])) for all pairs given above gives an (in most cases overdetermined) system of linear equations for the unknowncoefficients (Subscript[c, i]) . {(Subscript[y, 1]) , (Subscript[y, 2]) , ...,  (Subscript[y, s])}   is called the response vector. The values  (Subscript[ϕ, r])((Subscript[x, k]))  are the elements of the design matrix.  The latter is rectangular, in general. The number of columns is the number of functions, n. The number of rows is the number of points, s. Multiplying  the response vector with the pseudoinverse (cf. § 8.9 ) of the design matrix gives the unknown coefficients (Subscript[c, i]) .

Fit[data, funs, vars]  		finds a least square fit to a  list of  data  as a linear 
					combination of functions  funs  of variables in the list   vars

The function Fit[]  finds a least-squares fit of data to a linear combination of user selected basis functions.

Clear[data];
data = {{2,1}, {3,7}, {5,8}, {6,11}}; TableForm[data]

Clear[d,f,p]
d=ListPlot[data,PlotRange->{{0,6.1},{0,12}},Prolog->PointSize[0.015`]];
c=ListPlot[data,PlotRange->{{0,6.1},{0,12}},Joined->True];

Show[GraphicsRow[{d,c}]]

f = Fit[data, {1,x}, x]

pl = Plot[f, {x,0,6}, PlotRange -> {{0,6},{0, 12}}];

Clear[f]; f = Fit[data, {1,x,x^2}, x]

pq = Plot[f, {x,0,6}, PlotRange -> {{0,6},{0, 12}}];

ppl = Show[d,pl]; ppq = Show[d,pq];
 
Show[GraphicsRow[{ppl,ppq}]]

Clear[f]
f = Fit[data, {1, Sin[x]}, x]

p = Plot[f, {x,0,6}, PlotRange -> {{0,6},{0, 12}}];

Show[d,p, PlotRegion -> {{.01,.95},{.01,.95}}, ImageSize -> 200]

Comparison of the various curves produced, shows that their shapes depend very strongly on the functions one assumes. In particular, the extrapolation, i.e. the evaluation of values outside the interval, where the data were given, is extremely sensitive to these assumptions. Some general knowledge on the behaviour of the function one wants to obtain is necessary to get meaningful results.

FindFit[data, expr,par, var]      finds numerical values of the parameters par  so that
				   expr    gives a best fit to data.  var   lists the variable(s)

Clear[a,b,c,d,t, x, y]
FindFit[data,a + b t,{a,b},t]

data2 = Partition[{Range[Length[data]],Transpose[data]}//Flatten,3]

Map[Point,data2]

fxy = a + b  x + c y + d x y;
con =FindFit[data2,fxy,{a,b, c, d},{x, y}]

p1=Plot3D[Evaluate[fxy/.con],{x,0,8},{y,0,8},AxesLabel->{"x","y","fxy"},PlotLabel->Row[{"fxz  =  ",fxy}]];
p2 = Show[p1,Graphics3D[{PointSize[0.02`],Hue[0],Map[Point,data2]}], ImageSize -> 450]

12.2 Interpolation by Polynomials

InterpolatingPolynomial[data, var] 	gives a polynomial in the variable  var  which provides
				  			an exact fit to the data given in the list  data

Clear[c,f,g,p,x,y]; 
data = {{2,1}, {3,7}, {5,8}, {6,11}};TableForm[data]

poi = {AbsolutePointSize[8],Map[Point,data]};

f = InterpolatingPolynomial[data, x] //Expand

Plot[f,{x,0,6},PlotStyle -> Thickness[0.003],PlotRange -> {0,12}, Epilog -> poi]

For these four given points the interpolating polynomial gives a fit which is very satisfactory from an esthetic point 
of view. However, for larger sets of data points the curve representing the interpolating polynomials may oscillate 
very  strongly between the data points.

data1 = {{1,2},{2,1},{3,7},{4,6},{5,8},{5.5,7},{6,11}};
{PointSize[0.03],Map[Point,data1]};
d = Graphics[%,Axes -> True, PlotRange -> {{0,6.1},{0,12}}, AspectRatio -> 0.4];

g[x_] = Expand[InterpolatingPolynomial[data1,x]]

p=Plot[g[x],{x,0,6},PlotRange->{{0,6},{-8,12}}, Epilog -> d[[1]]];Show[GraphicsRow[{d,p}],ImageSize ->   500]

So one may prefer to use again least square fits employing a few powers only.

f = Fit[data1, {1,x}, x]

g = Fit[data1, {1,x,x^2}, x]

h = Fit[data1, {1,x,x^2,x^3}, x]

i = Fit[data1, {1,x,x^2,x^3,x^4}, x]

pf=Plot[f,{x,0,6},PlotLabel->"Linear Polynomial"];
pg=Plot[g,{x,0,6},PlotLabel->"Quadratic Polynomial"];
ph=Plot[h,{x,0,6},PlotLabel->"Cubic Polynom"];
pi=Plot[i,{x,0,6},PlotLabel->"Quartic Polynom"];
pdf=Show[pf,d];pdq=Show[pg,d];
pdh=Show[ph,d];pdi=Show[pi,d];
Show[GraphicsGrid[{{pdf,pdq},{pdh,pdi}}]]

12.3   Interpolation

??Interpolation

Function values and derivatives may be real or complex, or arbitrary symbolic expressions. The (Subscript[x, i]) must be real numbers.
Interpolation works by fitting polynomial curves between successive data points. The degree of these is specified by the
option  InterpolationOrder; whose default values is 3.  The precision of the InterpolatingFunction object is that of  the data.
Note: Lists comprising the data and the derivatives of the function(s) have been changed in going from Mathematica5.2
to Mathematica7.0 !

data = {-1,-0.5,0,1,3,6,9};Clear[lx,f]

f(x_)=((((E^(-x))-1))^2);

i1 = Interpolation[d1 = {data,f[data]}//Transpose];
epi =  {Epilog -> {Hue[1],PointSize[.025],Map[Point , d1]}};

p1=Plot[f[x],{x,-1,9},Evaluate[epi],PlotStyle->Hue[0.8`]];
p2=Plot[i1[t],{t,-1,9},Evaluate[epi]];
Show[GraphicsRow[{p1,p2}],ImageSize->500]

Using derivatives in interpolation

Comparison of the two drawings above shows that the interpolating function gives a curve lacking smoothness.
This can be improved by presribing also derivatives at the points.

d2 = Table[{{data[[k]]},f[data[[k]]],f'[data[[k]]]}//N,{k,Length[data]}]

Clear[i2,t]
i2[t_]=Interpolation[d2][t];

p3=Plot[i2[t],{t,-1,9},Evaluate[epi],PlotStyle->Hue[0.6], ImageSize ->   250]

Alternatively, one can use splines.

Show[p1,p2,p3,ImageSize->250]

FunctionInterpolation[expr, {x, xmin, xmax}] evaluates expr with x running from xmin to xmax and constructs an
	InterpolatingFunction       object which represents an approximate function corresponding to the result.

Simple Interpolation

In general, InterpolatingFunction objects are continuous but not smooth, as in the following example:

Needs["PlotLegends`"]

f[x_] = Exp[2x];

g = FunctionInterpolation[f[x],{x,-1,1}] ;

Plot[Evaluate[Table[D[g[x],{x,n}],{n,0,4}]],{x,-1,1},PlotRange->All,ImageSize->350,PlotStyle->Table[{Thick,Hue[(n/5)]},{n,0,4}],
PlotLegend->{"g(x)","g'(x)","g''(x)","(g^((3)))((x)","(g^((4)))(x)"},LegendPosition->{-1,0}]

It can be seen that these InterpolatingFunction objects are piecewise polynomial functions.

Smooth Interpolation

If the derivatives of the function are also available, one can use FunctionInterpolation[] as follows:

h = FunctionInterpolation[{f[x],f'[x],f''[x]},{x,-1,1}]

Plot[Evaluate[Table[(∂_({x,n}))h[x],{n,0,4}]],{x,-1,1},PlotRange->All,ImageSize->350,PlotStyle->Table[{Thick,Hue[(n/5)]},{n,0,4}],PlotLegend->{"h(x)","h'(x)","h''(x)","(h^((3)))(x)","(h^((4)))(x)"},LegendPosition->{-1,0}]

In general, the derivatives of an InterpolatingFunction object are continuous up to the minimum order
which is contained in its data.

Robert Knapp,  rknapp@wolrfram.com
TMJ, 8 (#3, 2002) 394/5.

12.4   Approximation by Spline Functions

Splines are polynomials interpolating between two data points. Depending on the degree of the spline also continuity
of higher derivatives is prescriped. In particular, cubic splines pass through all the points, they have continuous first
and second derivatives at all inner points. The second derivative is zero at the end points.
There are two packages for splines:  Graphics`Spline`  and   NumericalMath`SplineFit`. The first package calls and
uses the second one.

12.4.1  Drawing a spline curve through the data points.

<<Splines`;

data1={{1,2},{2,1},{3,7},{4,6},{5,8},{5.5`,7},{6,11}};
d = Graphics[{PointSize[0.01],Map[Point,data1]}, Axes -> True,
PlotRange->{{0,6},{0,15}}, AspectRatio ->   0.5]

Much more information and types of splines can be found  at  Help/  Find Selected Function /Splines .

sp=Show[Graphics[Spline[data1,Cubic]], AspectRatio -> 0.5];
dp=Show[d,sp];
Show[GraphicsRow[{sp,dp}],ImageSize->400]

12.4.2  Computing and using a spline

The data points need not be equidistant. For this demonstration we choose 4 data points:

np =  4 ;  (*  = number of data points *)

data=Table[{x+RandomReal[] 0.1`,RandomReal[]},{x,0,1,(1/(np-1))}]

dp=ListPlot[data,PlotRange->{{0,1.1`},{0,1.1`}},PlotStyle->PointSize[0.02`]];
sp=Graphics[Spline[data,Cubic]];
psd=Show[dp,sp,Axes->True];

Show[GraphicsRow[{dp,psd}],ImageSize->450]

A spline function object is obtained by the following command. It depends on a parameter, which runs
from 0 to np.  It may be used to plot a curve and to get values:

ip=SplineFit[data,Cubic];

pip=ParametricPlot[ip[u],{u,0,np-1},PlotRange->{0,1}, ImageSize ->   200, AspectRatio -> 0.6]

Gives the same curve as above at the right.  These are the coordinates of a point on the spline curve:

ip[2.23]

14.4.2.1   The splinefunction object

ip=SplineFit[data,Cubic];
fu = FullForm[ip]

The first argument characterizes the spline, the second one gives the range of the parameter, the third one 
the data points. The fourth argument gives for each of the np - 1 intervals bounded by two adjacent data points 
two lists; these are the coefficients  a0, a1, a2, a3  and b0, b1, b2, b3 
of cubic polynomials representing the spline curve in parametric form:

	{x[t] = a0  +  a1 t  +  a2 (t^2) + a3 (t^3), y[t] = b0  +  b1 t  +  b2 (t^2) + b3 (t^3)},    0 ≤  t  ≤  1.

ip[[4]]//Chop

tl = Table[t^n,{n,0,3}]

px = Table[ip[[4,n,1]].tl,{n,np-1}]//Chop

py = Table[ip[[4,n,2]].tl,{n,np-1}]//Chop

Properties of the splines

The above polynomials are used to show the properties of the splines at the data points.
They  are continuous at the inner points:

{px[[1]] , py[[1]] } /. t -> 0

Table[{(px[[k+1]] /. t -> 0) -( px[[k]] /. t -> 1), 
(py[[k+1]] /. t -> 0) - (py[[k]] /. t -> 1)}//Chop,{k,np - 2}]

{px[[-1]] , py[[-1]] }/. t -> 1

The first derivatives  are continuous at the inner points:

{D[px[[1]],t] , D[py[[1]],t] } /. t -> 0 //Chop

Table[{(D[px[[k+1]] ,t]/. t -> 0) -( D[px[[k]],t] /. t -> 1), 
(D[py[[k+1]],t] /. t -> 0) - (D[py[[k]],t] /. t -> 1)}//Chop,{k,np - 2}]

{D[px[[-1]],t] , D[py[[-1]],t] }/. t -> 1 //Chop

The second derivatives  are continuous at the inner points, zero at the end points:

{D[px[[1]],{t,2}] , D[py[[1]],{t,2}] } /. t -> 0 //Chop

Table[{(D[px[[k+1]] ,{t,2}]/. t -> 0) -( D[px[[k]],{t,2}] /. t -> 1), (D[py[[k+1]],{t,2}] /. t -> 0) - (D[py[[k]],{t,2}] /. t -> 1)}//Chop,{k,np - 2}]

{D[px[[-1]],{t,2}] , D[py[[-1]],{t,2}] }/. t -> 1//Chop

The first derivative of the spline curve

The first derivative is  dy/dx =  (OverDot[y])[t]/(OverDot[x])[t] ; the corresponding fractions are combined with x[t]
so that we get a parametric representation of the derivative of the spline curve:

fr = D[py,t]/D[px,t]

fd[u_] = Apply[Which,Flatten[Table[{k - 1 ≤  u ≤ k, 
{px[[k]], fr[[k]]} /. t -> u - k + 1}, {k,np-1} ]//ExpandAll//Together,1]]

pd=ParametricPlot[fd[u],{u,0,np-1},PlotStyle->Hue[0]];

Show[pip,pd,AxesLabel->{"x",""},PlotRange->{-3.5`,2.1`},Epilog->{Text["f(x)",{0.95`,0.6`}],Hue[0],Text["f'(x)",{0.5`,-0.9`}]}]

Second derivatives may be found from the formula: (d^2)y/d(x^2)  =  ((OverDot[x]) (Overscript[y, ..])  - (Overscript[x, ..])  (OverDot[y]))/((OverDot[x])^3).

12.5 Approximation by Chebyshev Expansions

A very economic method to represent functions is by expansions with respect to Chebyshev polynomials  Tn(x) (the German transliteration is:  Tschebischeff Polynome):

	  			Tn(x) : Tn(cos ϑ)  =  cos(n ϑ),   x = cos ϑ ,  -1 ≤ x ≤ 1;
	  			
	  			T0(x) = 1,   T1(x) = x,   T2(x)   =  (x^2)  - 1/2 , T3(x)   =  (x^3) - 3/ 4 x, ....

These fulfil the following orthonormality relation:

         					(∫_(-1)^(1  ))  dx Tn(x) Tk(x) / √(1 - x2)  =   δnk (1 +   δn0) π/2.

Tn(x)  is the polynomial of n-th degree with leading coefficient 1 whose maximum absolute value is smallest in
 the intervall   -1 ≤ x ≤ 1; i.e. it deviates least from zero.  In series expansions the shifted Chebyshev 
 polynomials are needed defined as:
							T*n(x)  =  Tn(2x - 1),  0 ≤ x ≤ 1;
with the following orthonormality properties:

(∫_0^(1  ))dx   T*n(x) T*k(x) / √(x - x2)  =   δnk (1 +   δn0) π/2.

The corresponding series exansion is:

					f(x)  =  Σ  cn  T*n(x),    0 ≤ x ≤ 1.

The advantage of the Chebyshev expansion is that it gives an approximation more uniform over the intervall, while Taylor's expansion is very accurate in the neighbourhood of the expansion point but becomes the less accurate the larger the distance from this point.  This is shown in a simple example  (due to Lanczos) where comparison with the exact result is easy.

f = 1/(1+x);

f6 = Series[f,{x,0,6}]//Normal

f7 = Series[f,{x,0,7}]//Normal

The expansion coefficients for the Chebyshev series are found as is usual for a complete orthonormal system. The necessary  integrations are done by  numeric integration. The factor  1/2  for   n = 0   is taken into account in a second step.

cn = Table[ NIntegrate[ChebyshevT[n,2x-1]/((1+x)Sqrt[x(1-x)]),{x,0,1}], 
		                    {n,0,6}] 2/Pi //N ;
cn[[1]] = cn[[1]]/2; cn

fc = Expand[Sum[ cn[[k]] ChebyshevT[k-1,2x-1],{k,Length[cn]}]];  fc

SetOptions[Plot,DisplayFunction->Identity];
pt=Plot[{f,fc,f6,f7},{x,0,1},AxesLabel->{"x","f, fc,f6,f7"},PlotStyle->{Dashing[{}],Dashing[{0.1`}],Dashing[{0.01`}],Dashing[{0.02`}]}];
pc=Plot[{fc-f},{x,0,1},PlotRange->(({-1,1})/(10^5)),AxesLabel->{"x","f - fc"}];
Show[GraphicsRow[{pt,pc}], ImageSize -> 500]

p6=Plot[{f-f6},{x,0,1},AxesLabel->{"x","f - f6"}];
p7=Plot[{f-f7},{x,0,1},AxesLabel->{"x","f - f7"}];
Show[GraphicsRow[{p6,p7},Spacings->{Scaled[0.5`],Scaled[0]}]]

f80=Normal[Series[f,{x,0,80}]];
SetOptions[Plot,AxesLabel->{"x","f - f80"}];
pt=Plot[{f-f80},{x,0,1},PlotRange->{0,-(1/(10^6))}];
pa=Plot[{f-f80},{x,0,1},PlotRange->All];
Show[GraphicsRow[{pt,pa},Spacings->{Scaled[0.5`],Scaled[0]}]]

Note that 80 (!) terms of the Taylors expansion give high precision up to x  ≤  .8, but thereafter the accuracy deterioriates rapidly.  -  In this simple example the expansion coefficients can be found by analytic integration;  in this way  one gets:

Clear[p]
fc = Sqrt[2] (1/2 + Sum[(-p)^n ChebyshevT[n,2x-1],{n,6}]); 
p = 3 - 2 Sqrt[2] ; N[p]
tc = N[fc]//Expand

Chop[fc - tc,10^-8]

dn = Table[
Integrate[ChebyshevT[n,2x-1]/((1+x) Sqrt[x(1-x)]),
	{x,0,1}], {n,0,6}] 2/Pi ;
dn[[1]] = dn[[1]]/2; Cancel[dn]
N[dn]

Prepend[Table[Sqrt[2] (-p)^k //Expand, {k,6}],(1/(Sqrt[2]))]

12.6 Nonlinear Fits

The function Fit[] contained in the kernel finds a least-squares fit of data to linear combination of given basis functions. Often a more complicated formula is needed, which is nonlinear in the parameters. The function  NonlinearModelFit[] can be used for that purpose.

??NonlinearModelFit

data={{0,1},{1,0},{3,2},{5,4}, {6,4},{7,5}};

nlm=NonlinearModelFit[data,Log[a +b x^2],{a,b},x]

Obtain the functional form:

Normal[nlm]

Evaluate the model at a point:

nlm[2.3]

Visualize the fitted function with the data:

Show[ListPlot[data],Plot[nlm[x],{x,0,7}],Frame->True]

Extract and plot the residuals:

nlm["FitResiduals"]

ListPlot[%,Filling->Axis]

{{x11,x12,...,y1}, {x21,x22,...,y2}, ...}  each sublist gives a data point.

The data set below consists of triples {x1,x2,y} describing a chemical reaction rate  y  as functions of two inputs x1 and x2.  The list gives 5 such sets of data.  The model is given by

y=((θ1 θ3 x1)/(x1 θ1+x2 θ2+1.))

The parameters θ1, θ2, θ3  must be determined by a least-squares fit from the following data set:

data =
	 {{1.,1.,.126},{1.,2.,.076},{2.,1.,.219},{2.,2.,.126},{.1, 0.,.186}};

nlm = NonlinearModelFit[data,(θ1*θ3*x1)/(x1*θ1 + x2*θ2 + 1),{θ1, θ2, θ3},{x1,x2}]

Normal[nlm]

12.7 Exercises

12.1. 	Find a curve composed of cubic splines for the set of data named data in   section 12.1 .
     
12.2. 	Generate a list of 5 complex random numbers located in the square, whose corners are at the origin and  at  5 + 5i. 
	Tranform each of this complex numbers in a pair of real numbers; each represents a point in the real plane. 
	Fit a polynomial of first, second, third degree and a curve of cubic splines to these points. 
	Draw each of these curves together with the data points.

12.3  	Find a Chebyshev expansion comprising 7 terms for  sin(x) in  the intervall 0  ≤  x   ≤  π/2; compare its accuracy to 
	that of a Taylor's expansions performed around  x = 0 and  x = π/4.

12.4	Fit the data given below to the implicit function 1/x = exp( b y/ x ). 	Determine b such that it gives the best fit.

data = {{.5, -.693147}, {1.,0.}, {2.,2.7759}, {3.,6.592}};

Hint: Consider a function   f = 1/x - exp( b y/ x ) = 0. Augment data to new list data1 for {x,y,f} as input for  
	NonLinearModelFit[].

12.5  	Plot the surface corresponding to the final model function given at the end of sect.12.6
	together with points representing the input data.
	Compute the residuals and show them together with the curcumference of a square at height  h.

12.6	A set of points is given in 3 dimensions. Find the best fit plane determined by these points.
	The equation of this plane is: a x + b y + c z = d; find values for a, b, c, d.