Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Least Squares & Gram-Schmidt to Construct Orthogonal Polynomials

Needed: Norm, Guassian elimination module, an inverse matrix module

norm2[g_]:=((((∑_(i=1)^(Length[g]))((g[[i]])^2)))^(1/2));

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Return[A];  ] ;

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];A3=Take[A2,{1,n},{n+1,2n}];Return[A3]];

Linear Least Squares: Normal Equations

Data Set

Here is the data set:

data={{4.0,102.56},{4.2,113.18},{4.5,130.11},{4.7,142.05},{5.1,167.53},{5.5,195.14},{5.9,224.87},{6.3,256.73},{6.8,299.50},{7.1,326.72}} ;

n=Length[data];

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]}]

Matrix Model y=(Subscript[c, 1])+(Subscript[c, 2])x

The model is y=(Subscript[c, 1])+(Subscript[c, 2])x

Clear[x,y,f1,f2,A,a,B,b,AT,at]
(*Define the linear model that will be fit to the data*)
m=2;
f1[x_]=1;
f2[x_]=x;

(*Construct the matrix A*)
A=Array[a,{n,m}];
Do[
a[j,1]=f1[data[[j,1]]];
a[j,2]=f2[data[[j,1]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{n}];
Do[
b[j]=data[[j,2]];
,{j,1,n}];
MatrixForm[B]

(*Construct the Transpose matrix*)
AT=Array[at,{m,n}];
Do[
at[j,i]=a[i,j];
,{j,1,m}
,{i,1,n,1}];
MatrixForm[AT]

The normal equations

ATA=AT.A
ATB=AT.B

Put them into a matrix that we can send to the Gaussian Elimination Module to solve.

matrix1=SparseArray[{{0}},{m,m}];
Table[matrix1[[i]]=Append[ATA[[i]],ATB[[i]]],{i,1,m}];
MatrixForm[matrix1]

Use the Gaussian Elimination module to solve the system:

sol=Array[f,m];
Do[sol[[i]]=GaussElim3[matrix1][[i,m+1]],{i,1,m}];
sol

Define the value of the parameters:

c1=sol[[1]]
c2=sol[[2]]
y[x_]=c1*f1[x]+c2 *f2[x]

Plot the least squares fit!

model = Plot[y[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{{4.0,7.1},{90,400}}]

The measure of the error in the fit can be the squared error:

SE=Sum[(data[[i,2]]-y[data[[i,1]]])^2,{i,1,n}]
(Sqrt[SE])

which we can also compute as follows:

B-A.sol
norm2[B-A.sol]

or the root mean square error:

RMSE=(Sqrt[SE/Length[data]])

Matrix Model y=(Subscript[c, 1])+(Subscript[c, 2])x+(Subscript[c, 3])(x^2)

The model is y=(Subscript[c, 1])+(Subscript[c, 2])x+(Subscript[c, 3])(x^2)

Clear[x,y,f1,f2,f3,A,a,B,b,AT,at]
(*Define the linear model that will be fit to the data*)
m=3; (*This is the number of parameters*)
f1[x_]=1;
f2[x_]=x;
f3[x_]=(x^2);

(*Construct the matrix A*)
A=Array[a,{n,m}];
Do[
a[j,1]=f1[data[[j,1]]];
a[j,2]=f2[data[[j,1]]];
a[j,3]=f3[data[[j,1]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{n}];
Do[
b[j]=data[[j,2]];
,{j,1,n}];
MatrixForm[B]

(*Construct the Transpose matrix*)
AT=Array[at,{m,n}];
Do[
at[j,i]=a[i,j];
,{j,1,m}
,{i,1,n,1}];
MatrixForm[AT]

The normal equations

ATA=AT.A
ATB=AT.B

Put them into a matrix that we can send to the Gaussian Elimination Module to solve.

matrix1=SparseArray[{{0}},{m,m}];
Table[matrix1[[i]]=Append[ATA[[i]],ATB[[i]]],{i,1,m}];
MatrixForm[matrix1]

Use the Gaussian Elimination module to solve the system:

sol=Array[f,m];
Do[sol[[i]]=GaussElim3[matrix1][[i,m+1]],{i,1,m}];
sol

Define the value of the parameters:

c1=sol[[1]]
c2=sol[[2]]
c3=sol[[3]]
y[x_]=c1*f1[x]+c2 *f2[x]+c3*f3[x]

Plot the least squares fit!

model = Plot[y[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{{4.0,7.1},{90,400}}]

The measure of the error in the fit can be the squared error:

SE=Sum[(data[[i,2]]-y[data[[i,1]]])^2,{i,1,n}]

or the root mean square error:

RMSE=(Sqrt[SE/Length[data]])

Matrix Model y=(Subscript[c, 1])(x^(Subscript[c, 2]))

In this case, the setup is identical to that of the polynomial fits, except we have to preprocess the fit function f(x)=b (x^a) to create a linear function. We can write

y=(Subscript[c, 1]) (x^(Subscript[c, 2]))
ln y = ln (Subscript[c, 1]) + (Subscript[c, 2]) ln x
(OverTilde[y])=ln (Subscript[c, 1])+(Subscript[c, 2]) (OverTilde[x])

So we have to convert the data by taking logarithms.

Clear[x,y,f1,f2,A,a,B,b,AT,at]
(*Define the linear model that will be fit to the data*)
m=2;
f1[x_]=1;
f2[x_]=x;

(*Construct the matrix A*)
A=Array[a,{n,m}];
Do[
a[j,1]=f1[Log[data[[j,1]]]];
a[j,2]=f2[Log[data[[j,1]]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{n}];
Do[
b[j]=Log[data[[j,2]]];
,{j,1,n}];
MatrixForm[B]

(*Construct the Transpose matrix*)
AT=Array[at,{m,n}];
Do[
at[j,i]=a[i,j];
,{j,1,m}
,{i,1,n,1}];
MatrixForm[AT]

The normal equations

ATA=AT.A
ATB=AT.B

Put them into a matrix that we can send to the Gaussian Elimination Module to solve.

matrix1=SparseArray[{{0}},{m,m}];
Table[matrix1[[i]]=Append[ATA[[i]],ATB[[i]]],{i,1,m}];
MatrixForm[matrix1]

Use the Gaussian Elimination module to solve the system:

sol=Array[f,m];
Do[sol[[i]]=GaussElim3[matrix1][[i,m+1]],{i,1,m}];
sol

Define the value of the parameters:

c1=Exp[sol[[1]]]
c2=sol[[2]]
modelfcn[x_]=c1 (x^c2)

Plot the least squares fit!

model = Plot[modelfcn[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{{4.0,7.1},{90,400}}]

The measure of the error in the fit can be the squared error:

SE=Sum[(data[[i,2]]-modelfcn[data[[i,1]]])^2,{i,1,Length[data]}]

or the root mean square error:

RMSE=(Sqrt[SE/Length[data]])

Modified Gram-Schmidt Orthogonalization: Reduced and Full QR Decomposition

Given m linearly independent vectors, we can construct mutually independent perpendicular vectors spanning the same space as the original m vectors.

set1={{1,-4},{2,3},{2,2}}
MatrixForm[set1]

To access the column vector of the matrix we can use the following command:

Transpose[set1][[1]]

Here is the modified Gram-Schmidt algorithm to do the QR decomposition:

Clear[n,m,Q,R,r,y]
n=Length[set1[[1]]];
m=Length[set1];
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
y=Transpose[set1][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].y;
y=y-r[i,j]q[i];
}];
r[j,j]=norm2[y];
q[j]=y/r[j,j];
Q=Append[Q,q[j]];
}];

MatrixForm[Simplify[Q]]
MatrixForm[Simplify[R]]

Here is the built-in Mathematica command to do the QR decomposition:

res=QRDecomposition[set1];
MatrixForm[res[[1]]]
MatrixForm[res[[2]]]//Simplify

Notice that our result matches what Mathematica returns, but to match the notation of column vectors in the text we should take the transpose of the matrix Q at the end. This is annoying, but for the moment I think it is more important to match the notation of the text.

Clear[n,m,Q,R,r,y]
n=Length[set1[[1]]];
m=Length[set1];
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
y=Transpose[set1][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].y;
y=y-r[i,j]q[i];
}];
r[j,j]=norm2[y];
q[j]=y/r[j,j];
Q=Append[Q,q[j]];
}];
Q=Transpose[Q]; (*added to match notation of the text*)

MatrixForm[Simplify[Q]]
MatrixForm[Simplify[R]]

Now this matches the text. And we can see that we did the factoring correctly, since A=Q R (with our fudge to transform Q at the end):

set1//MatrixForm
Q.R//MatrixForm

What we have so far is the reduced QR factorization. Note that until we go to the full QR decomposition, we do not have the important property we need for least squares, that (Q^T)=(Q^(-1)) since the inverse is only defined on a square matrix.

Transpose[Q].Q//MatrixForm

Full QR decomposition pads everything out to be square by adding an extra vector (any vector that is not linearly dependent to the others will do), and that is easiest to do when you construct the matrix.

set1={{1,-4,1},{2,3,0},{2,2,0}}
MatrixForm[set1]

Clear[n,m,Q,R,r,y]
n=Length[set1[[1]]];
m=Length[set1];
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
y=Transpose[set1][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].y;
y=y-r[i,j]q[i];
}];
r[j,j]=norm2[y];
q[j]=y/r[j,j];
Q=Append[Q,q[j]];
}];
Q=Transpose[Q]; (*added to match notation of the text*)

MatrixForm[Simplify[Q]]
MatrixForm[Simplify[R]]

Now we have the important property we need for least squares, that (Q^T)=(Q^(-1)):

set1//MatrixForm
Q.R//MatrixForm
Transpose[Q].Q//MatrixForm

Least Squares by Full QR Decomposition

Rather than using Gaussian elimination, we can use QR Decomposition to do linear least squares.

Here is the data set:

data={{4.0,102.56},{4.2,113.18},{4.5,130.11},{4.7,142.05},{5.1,167.53},{5.5,195.14},{5.9,224.87},{6.3,256.73},{6.8,299.50},{7.1,326.72}} ;

numpts=Length[data];

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]}]

We need to make the original matrix A be an n×n matrix, so we pad it out. Note that n is the number of data points, and m is the number of parameters (Subscript[c, i]) in our model. In this example, only the first two columns are intialized from the data set, the rest are padding. The padding is built from the identity matrix to ensure linear independence.

Clear[x,y,f1,f2,A,a,B,b,c,AT,at,n,m]
(*Define the linear model that will be fit to the data*)
m=2;
n=Length[data]
f1[x_]=1.0;
f2[x_]=x;

(*Construct the matrix A*)
A=IdentityMatrix[{n,n}];
Do[
A[[j,1]]=f1[data[[j,1]]];
A[[j,2]]=f2[data[[j,1]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{Length[data]}];
Do[
b[j]=data[[j,2]];
,{j,1,n}];
MatrixForm[B]

Here is the modified Gram-Schmidt algorithm to do the full QR decomposition.We will always want to be doing full QR decomposition on an n×n matrix (the total number of data points) when we are using it for linear least squares. I’ve also used a different y since that variable is often used in models.

Clear[Q,R,r,yQR]
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
yQR=Transpose[A][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].yQR;
yQR=yQR-r[i,j]q[i];
}];
r[j,j]=norm2[yQR];
q[j]=yQR/r[j,j];
Q=Append[Q,q[j]];
}];
Q=Transpose[Q]; (*added to match notation of the text*)

Here is what the full QR decomposition gives us:

R//MatrixForm
Transpose[Q].B//MatrixForm

Take only the bit that we need:

Rhat=Take[R,m,m];
MatrixForm[Rhat]
dhat=Take[Transpose[Q].B,m];
MatrixForm[dhat]

Solve the system (OverHat[R ])c = (OverHat[d]) by backsubstitution (since (OverHat[R]) is upper triangular).

Do[c[i]=((dhat[[i]]-(∑_(j=i+1)^m)c[j]Rhat[[i,j]])/(Rhat[[i,i]])),{i,m,1,-1}];

Here is our model, which is what we found before when we solved the normal equations.

modelfcn[x_]=c[1]*f1[x]+c[2]*f2[x]

model = Plot[modelfcn[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{{4.0,7.1},{90,400}}]

SE=Sum[(data[[i,2]]-modelfcn[data[[i,1]]])^2,{i,1,n}]

Linear Least Squares: More Complicated Model

Data Set

Here is a data set:

data={{0.1063477898034214386,2.0251095726333425695},{0.6057905198926769597,2.0128198521606748272},{1.0948272193491530579,1.9769560032250161896},{1.6008449324101925448,1.9337635859029760343},{2.0979201216309717459,1.8901566366681423759},{2.5942365195810825895,1.8343195196396289148},{3.0984570129867989863,1.769111429105063138},{3.6095434352348156892,1.7467646098369875585},{4.1085055040235767265,1.7175035620349641263},{4.5958573947359582235,1.7018748416582440974},{5.1016094896258827873,1.7096615005319446556},{5.5961285476006921158,1.746444294921263074},{6.0978116374276014611,1.7730557440898282451},{6.5965070299599055659,1.7883293637461021852},{7.1003754833408025959,1.7919725586821759755},{7.5994513028639287811,1.8164984407820472523},{8.0902367092559668293,1.8008386224407406662},{8.6057960805376136904,1.782746607219734436},{9.0978255277999213746,1.7811649629750643253},{9.5991768976590634566,1.7737489372743693306}};

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]}]

You have an idea that the function that models this looks something like (y^2)=(Subscript[c, 1])sin(x)/x +(Subscript[c, 2]), but you are not sure what values of the parameters c to choose. Let’s use QR to determine the (Subscript[c, i])

QR

The first thing is to decide how you want to write your model. If we choose:

(y^2)=(Subscript[c, 1])((sin(x))/x)+(Subscript[c, 2])
then
(Subscript[f, 1])(x) = ((sin(x))/x)
(Subscript[f, 2])(x)=1.0

and the vector B will have elements ((((Subscript[y, i])))^2)

Clear[x,y,f1,f2,A,a,B,b,c,n,m]
(*Define the linear model that will be fit to the data*)
m=2;
n=Length[data];
f1[x_]=Sin[x]/x;
f2[x_]=1.0;

(*Construct the matrix A*)
A=IdentityMatrix[{n,n}];
Do[
A[[j,1]]=f1[data[[j,1]]];
A[[j,2]]=f2[data[[j,1]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{n}];
Do[
b[j]=((data[[j,2]])^2);
,{j,1,n}];
MatrixForm[B]

Here is the modified Gram-Schmidt algorithm to do the full QR decomposition.

Clear[Q,R,r,yQR]
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
yQR=Transpose[A][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].yQR;
yQR=yQR-r[i,j]q[i];
}];
r[j,j]=norm2[yQR];
q[j]=yQR/r[j,j];
Q=Append[Q,q[j]];
}];
Q=Transpose[Q]; (*added to match notation of the text*)

Here is what the full QR decomposition gives us:

R//MatrixForm
Transpose[Q].B//MatrixForm

Take only the bit that we need:

Rhat=Take[R,m,m];
MatrixForm[Rhat]
dhat=Take[Transpose[Q].B,m];
MatrixForm[dhat]

Solve the system (OverHat[R ])c = (OverHat[d]) by backsubstitution (since (OverHat[R]) is upper triangular).

Do[c[i]=((dhat[[i]]-(∑_(j=i+1)^m)c[j]Rhat[[i,j]])/(Rhat[[i,i]])),{i,m,1,-1}];

Here is our model

modelfcn[x_]=(Sqrt[c[1]*f1[x]+c[2]*f2[x]])

model = Plot[modelfcn[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{All,All}]

And the error in the model is

SE=Sum[(data[[i,2]]-modelfcn[data[[i,1]]])^2,{i,1,n}]

QR Same Data, Different Model

If we choose:

(Subscript[c, 1])x (y^2)+(Subscript[c, 2])x=sin(x)
then
(Subscript[f, 1])(x,y) = x (y^2)
(Subscript[f, 2])(x)=x

and the vector B will have elements sin((Subscript[y, i]))

Clear[x,y,f1,f2,A,a,B,b,c,n,m]
(*Define the linear model that will be fit to the data*)
m=2;
n=Length[data];
f1[x_,y_]=x*(y^2);
f2[x_]=x;

(*Construct the matrix A*)
A=IdentityMatrix[{n,n}];
Do[
A[[j,1]]=f1[data[[j,1]],data[[j,2]]];
A[[j,2]]=f2[data[[j,1]]];
,{j,1,n}];
MatrixForm[A]

(*Construct the vector B*)
B=Array[b,{n}];
Do[
b[j]=Sin[data[[j,1]]];
,{j,1,n}];
MatrixForm[B]

Here is the modified Gram-Schmidt algorithm to do the full QR decomposition.

Clear[Q,R,r,yQR]
Q={};
R=Array[r,{n,n}];
Do[r[i,j]=0,{i,1,n},{j,1,i}]

For[j=1,j≤n,j=j+1,
{
yQR=Transpose[A][[j]];
For[i=1,i≤j-1,i=i+1,
{
r[i,j]=q[i].yQR;
yQR=yQR-r[i,j]q[i];
}];
r[j,j]=norm2[yQR];
q[j]=yQR/r[j,j];
Q=Append[Q,q[j]];
}];
Q=Transpose[Q]; (*added to match notation of the text*)

Here is what the full QR decomposition gives us:

R//MatrixForm
Transpose[Q].B//MatrixForm

Take only the bit that we need:

Rhat=Take[R,m,m];
MatrixForm[Rhat]
dhat=Take[Transpose[Q].B,m];
MatrixForm[dhat]

Solve the system (OverHat[R ])c = (OverHat[d]) by backsubstitution (since (OverHat[R]) is upper triangular).

Do[c[i]=((dhat[[i]]-(∑_(j=i+1)^m)c[j]Rhat[[i,j]])/(Rhat[[i,i]])),{i,m,1,-1}];

c[1]

c[2]

Here is our model

modelfcn[x_,y_]=c[1]*f1[x,y]+c[2]*f2[x]-Sin[x]

model = ContourPlot[modelfcn[x,y]==0,{x,0,10},{y,1.6,2},ContourStyle->{Thick,Darker[Green]}];
compplot2 = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{All,All}]

Show[compplot,compplot2]

Notice our two models are slightly different.

Nonlinear Least Squares: Newton’s Method

Data Set

Here is the data set:

data={{4.0,102.56},{4.2,113.18},{4.5,130.11},{4.7,142.05},{5.1,167.53},{5.5,195.14},{5.9,224.87},{6.3,256.73},{6.8,299.50},{7.1,326.72}} ;

n=Length[data];

dataplot=ListPlot[data,PlotStyle->{Red,PointSize[0.02]}]

Matrix Model y=(Subscript[c, 1])(x^(Subscript[c, 2]))

Let’s construct the system of equations to solve based on our model.

f(x)=(Subscript[c, 1]) (x^(Subscript[c, 2]))
(∂/(∂(Subscript[c, 1])))f(x)= (x^(Subscript[c, 2]))
(∂/(∂(Subscript[c, 2])))f(x)=ln(x) (Subscript[c, 1]) (x^(Subscript[c, 2]))

Clear[f,x,y,xp,yp,c1,c2,n,m]
n=Length[data];
m=2;
f[x_,c1_,c2_]=c1 x^c2;

D[f[x,c1,c2],c1]
D[f[x,c1,c2],c2]

Define the data points

Do[
xp[i]=data[[i,1]];
yp[i]=data[[i,2]];
,{i,1,n}]

eq1=Simplify[Sum[(yp[i]-f[xp[i],c1,c2])D[f[xp[i],c1,c2],c1],{i,1,n}]]
eq2=Simplify[Sum[(yp[i]-f[xp[i],c1,c2])D[f[xp[i],c1,c2],c2],{i,1,n}]]

We need to solve the equations using Newton’s method. Let’s make sure the Jacobian is defined before we begin.

f1[c1_,c2_]={eq1,eq2};
jacobian[A0_,vars_]:=Outer[D,A0,vars]
jac[c1_,c2_]=jacobian[f1[c1,c2],{c1,c2}];

Here we go! You need to choose a good starting value, so if you have to sample the solution space and find a place to start.

Clear[p0,p1,tol, max, i,p,q]

f1[c 1_,c2_]={eq1,eq2};

p0={5.0,3.0};
tol=10.0^(-10);
max=20;

data={};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

y= MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p = p0 -y;

data=Append[data,{k,p,norm2[y]}];
If[norm2[y]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","Norm p"}}]

So we have our solutions. The model is

modelfcn[x_]=6.2422646831939028189*(x^2.0192443927498047174)

model = Plot[modelfcn[x],{x,0,10},PlotStyle->{Thick,Blue}];
compplot = Show[model,dataplot,Frame->True,FrameLabel->{"t","y"},PlotRange->{{4.0,7.1},{90,400}}]

The measure of the error in the fit can be the squared error:

SE=Sum[(yp[i]-modelfcn[xp[i]])^2,{i,1,n}]

or the root mean square error:

RMSE=(Sqrt[SE/Length[data]])

Aside: Least Squares: Continuous

Constructing Legendre Polynomials: Gram-Schmidt Process

Given an interval and a weight function, we are going to construct an orthogonal polynomial using the Gram-Schmidt process. This gives us another way of obtaining orthogonal polynomial that does not depend on the series solution to a differential equation!

Here, the weight function is w(x)=1and the interval is [-1,1]. We already know the orthogonal polynomials associated with this weight and interval are the Legendre polynomials, so that's what we should find.

Clear[ϕ,w,B,A,a,b,num,k,n]

a=-1;
b=1;
w[x_]=1;

(Subscript[ϕ, 0])[x_]=1;
num=4;

(Subscript[B, k_]):=(Subscript[B, k])=(((∫_a^b)x*w[x]*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x)/((∫_a^b)w[x]*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x))

(Subscript[A, k_]):=(Subscript[A, k])=(((∫_a^b)x*w[x]*(Subscript[ϕ, k-1])[x]*(Subscript[ϕ, k-2])[x]\[DifferentialD]x)/((∫_a^b)w[x]*(((Subscript[ϕ, k-2])[x])^2)\[DifferentialD]x))

(Subscript[ϕ, 1])[x_]=x-(Subscript[B, 1]);

Do[(Subscript[ϕ, k])[x_]=Expand[(x-(Subscript[B, k]))(Subscript[ϕ, k-1])[x]-(Subscript[A, k])(Subscript[ϕ, k-2])[x]],{k,2,num}]
TableForm[Table[{n,(Subscript[ϕ, n])[x]},{n,0,num}]]

TableForm[Table[{k,LegendreP[k,x]},{k,0,num}]]

We are missing the standardization (Subscript[P, n])(x)=1, but other than that, what we have found is the Legendre Polynomials!

Curve Fitting

Here is a function we wish to approximate using thge Legendre Polynomials:

f[x_]=2/(x^2+x+2)
fcnplot=Plot[f[x],{x,-2.,2},Frame->True,PlotStyle->{Red,Thick}]

Clear[approx,a]
(Subscript[β, k_]):=(Subscript[β, k])=(2/(2 k+1))
(Subscript[a, k_]):=(Subscript[a, k])=((NIntegrate[LegendreP[k,x] f[x],{x,-1,1}])/(Subscript[β, k]))
approx[n_,x_]:=approx[n,x]=(∑_(i=0)^n)(Subscript[a, i]) LegendreP[i,x]
Table[{n,N[approx[n,x]]},{n,0,3}]

plot[n_]:=Plot[{f[x],approx[n,x]},{x,-2,2},PlotStyle->{{Red,Thick},{Blue,Thick}},PlotRange->{{-1.4,1.4},{0.25,1.25}}]

Show[GraphicsRow[{plot[1],plot[2],plot[3]},DisplayFunction->$DisplayFunction,ImageSize->{650,Automatic}]]

Here is an estimate of the error:

TableForm[Table[{n,NIntegrate[(f[x]-approx[n,x])^2,{x,-1,1}]},{n,0,3}]]