(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = leftheader, inactive, L2,  12, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; 
	fontset = leftfooter, inactive, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	paletteColors = 128; currentKernel; 
]
:[font = title; inactive; Cclosed; preserveAspect; startGroup]
Approximation of Functions-part iii
:[font = subtitle; inactive; preserveAspect]
An interpretation based upon Theodore J. Rivilin's "An Introduction to the Approximation of Functions"
:[font = subsubtitle; inactive; preserveAspect; endGroup]
M.A. Lachance
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Classical orthogonal polynomials
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Chebyshev polynomials
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Chebyshev polnomials of the FIRST kind
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[f,g,w,a,b,alpha,beta];
w[x_]:= 1/Sqrt[1-x^2];
a = -1;
b=1;

f[0,x_]=1;

alpha[0]= Integrate[x*w[x],{x,a,b}]/
          Integrate[ w[x],{x,a,b}];
f[1,x_]=x-alpha[0];

For[i=2,i<=5,i++,
    alpha[i-1]= 
      Integrate[x*f[i-1,x]*f[i-1,x]*w[x],{x,a,b}]/
        Integrate[ f[i-1,x]*f[i-1,x]*w[x],{x,a,b}];
    beta[i-1]= 
      Integrate[x*f[i-1,x]*f[i-2,x]*w[x],{x,a,b}]/
        Integrate[ f[i-2,x]*f[i-2,x]*w[x],{x,a,b}];
    f[i,x_]=Expand[(x-alpha[i-1])f[i-1,x]-beta[i-1] f[i-2,x]];
    ];
;[s]
3:0,0;4,1;19,0;521,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect]
MatrixForm[Table[Expand[f[n,x]/f[n,1]],{n,0,5}]]
MatrixForm[Table[Expand[ChebyshevT[n,x]/
ChebyshevT[n,1]],{n,0,5}]]
:[font = input; wordwrap; preserveAspect; endGroup]
Plot[{
ChebyshevT[0,x],
ChebyshevT[1,x],
ChebyshevT[2,x],
ChebyshevT[3,x],
ChebyshevT[4,x],
ChebyshevT[5,x],
ChebyshevT[6,x],
ChebyshevT[7,x],
ChebyshevT[8,x],
ChebyshevT[9,x]},{x,-1,1},
AspectRatio->Automatic]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Chebyshev polnomials of the SECOND kind
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[f,g,w,a,b,alpha,beta];
w[x_]:= Sqrt[1-x^2];
a = -1;
b=1;

f[0,x_]=1;

alpha[0]= Integrate[x*w[x],{x,a,b}]/
          Integrate[ w[x],{x,a,b}];
f[1,x_]=x-alpha[0];

For[i=2,i<=5,i++,
    alpha[i-1]= 
      Integrate[x*f[i-1,x]*f[i-1,x]*w[x],{x,a,b}]/
        Integrate[ f[i-1,x]*f[i-1,x]*w[x],{x,a,b}];
    beta[i-1]= 
      Integrate[x*f[i-1,x]*f[i-2,x]*w[x],{x,a,b}]/
        Integrate[ f[i-2,x]*f[i-2,x]*w[x],{x,a,b}];
    f[i,x_]=Expand[(x-alpha[i-1])f[i-1,x]-beta[i-1] f[i-2,x]];
    ];
;[s]
3:0,0;4,1;19,0;519,-1;
2:2,7,10,Courier,1,12,0,0,0;1,7,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect]
MatrixForm[Table[Expand[f[n,x]/f[n,1]],{n,0,5}]]
MatrixForm[Table[Expand[ChebyshevU[n,x]/
ChebyshevU[n,1]],{n,0,5}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Notice "arcsine envelope"
:[font = input; wordwrap; preserveAspect; endGroup]
Plot[{
ChebyshevU[0,x],
ChebyshevU[1,x],
ChebyshevU[2,x],
ChebyshevU[3,x],
ChebyshevU[4,x],
ChebyshevU[5,x],
ChebyshevU[6,x],
ChebyshevU[7,x],
ChebyshevU[8,x],
ChebyshevU[9,x]},{x,-1,1},
AspectRatio->Automatic]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Multiplying by SQRT[weight] produces nearly equioscillatory behavior
:[font = input; wordwrap; preserveAspect; endGroup; endGroup]
Plot[{
Sqrt[w[x]]ChebyshevU[0,x],
Sqrt[w[x]]ChebyshevU[1,x],
Sqrt[w[x]]ChebyshevU[2,x],
Sqrt[w[x]]ChebyshevU[3,x],
Sqrt[w[x]]ChebyshevU[4,x],
Sqrt[w[x]]ChebyshevU[5,x],
Sqrt[w[x]]ChebyshevU[6,x],
Sqrt[w[x]]ChebyshevU[7,x],
Sqrt[w[x]]ChebyshevU[8,x],
Sqrt[w[x]]ChebyshevU[9,x]},{x,-1,1},
AspectRatio->Automatic]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
The first and second kinds are the linearly indpendent solutions of a second order ODE
:[font = input; wordwrap; preserveAspect; startGroup]
n=5;

Clear[f,g];
f[x_] := ChebyshevT[n,x];
g[x_] := Sqrt[1-x^2] ChebyshevU[n-1,x];
Simplify[Expand[(1-x^2) f''[x]- x f'[x] + n^2 f[x]]]
Simplify[Expand[(1-x^2) g''[x]- x g'[x] + n^2 g[x]]]
:[font = output; output; inactive; preserveAspect]
0
;[o]
0
:[font = output; output; inactive; preserveAspect; endGroup; endGroup; endGroup]
0
;[o]
0
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Legendre functions
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Legendre polnomials of the FIRST kind
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[f,g,w,a,b,alpha,beta];
w[x_]:= 1;
a = -1;
b=1;

f[0,x_]=1;

alpha[0]= Integrate[x*w[x],{x,a,b}]/
          Integrate[ w[x],{x,a,b}];
f[1,x_]=x-alpha[0];

For[i=2,i<=5,i++,
    alpha[i-1]= 
      Integrate[x*f[i-1,x]*f[i-1,x]*w[x],{x,a,b}]/
        Integrate[ f[i-1,x]*f[i-1,x]*w[x],{x,a,b}];
    beta[i-1]= 
      Integrate[x*f[i-1,x]*f[i-2,x]*w[x],{x,a,b}]/
        Integrate[ f[i-2,x]*f[i-2,x]*w[x],{x,a,b}];
    f[i,x_]=Expand[(x-alpha[i-1])f[i-1,x]-beta[i-1] f[i-2,x]];
    ];
;[s]
3:0,0;4,1;19,0;509,-1;
2:2,7,10,Courier,1,12,0,0,0;1,7,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect]
MatrixForm[Table[Expand[f[n,x]/f[n,1]],{n,0,5}]]
MatrixForm[Table[Expand[LegendreP[n,x]/
LegendreP[n,1]],{n,0,5}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Nearly equioscillatory behavior
:[font = input; wordwrap; preserveAspect; endGroup; endGroup]
Plot[{
LegendreP[0,x],
LegendreP[1,x],
LegendreP[2,x],
LegendreP[3,x],
LegendreP[4,x],
LegendreP[5,x],
LegendreP[6,x],
LegendreP[7,x],
LegendreP[8,x],
LegendreP[9,x]},{x,-1,1},
AspectRatio->Automatic]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Legendre fucntions of the SECOND kind
:[font = input; preserveAspect; endGroup]
MatrixForm[Table[LegendreQ[i1,x],{i1,0,4}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
The first and second kinds are the linearly indpendent solutions of a second order ODE
:[font = input; wordwrap; preserveAspect; endGroup; endGroup]
n=5;

Clear[f,g];
f[x_]:= LegendreP[n,x];
g[x_]:= LegendreQ[n,x];

Simplify[Expand[
(1-x^2) f''[x]- 2 x f'[x]+n(n+1)f[x]]]

Simplify[Expand[
(1-x^2) g''[x]-2 x g'[x]+n(n+1)g[x]]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Laguerre polnomials
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[f,g,w,a,b,alpha,beta];
w[q_,x_]:= x^q Exp[-x];
q = 1;
a = 0;
b = Infinity;

f[0,x_]=1;

alpha[0]= Integrate[x*w[q,x],{x,a,b}]/
          Integrate[ w[q,x],{x,a,b}];
f[1,x_]=x-alpha[0];

For[i=2,i<=5,i++,
    alpha[i-1]= 
      Integrate[Expand[x*f[i-1,x]*f[i-1,x]]*w[q,x],{x,a,b}]/
        Integrate[ Expand[f[i-1,x]*f[i-1,x]]*w[q,x],{x,a,b}];
    beta[i-1]= 
      Integrate[Expand[x*f[i-1,x]*f[i-2,x]]*w[q,x],{x,a,b}]/
        Integrate[Expand[ f[i-2,x]*f[i-2,x]]*w[q,x],{x,a,b}];
    f[i,x_]=Expand[(x-alpha[i-1])f[i-1,x]-beta[i-1] f[i-2,x]];
    ];
;[s]
3:0,0;4,1;19,0;581,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect]
MatrixForm[Table[Expand[f[n,x]/f[n,1]],{n,0,5}]]
MatrixForm[Table[Expand[LaguerreL[n,q,x]/
LaguerreL[n,q,1]],{n,0,5}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Nearly equioscillatory behavior
:[font = input; wordwrap; preserveAspect; endGroup; endGroup]
q = 20;
n = 50;
Plot[{
  Sqrt[w[q, x]] LaguerreL[n, q, x]}, 
  {x, .01, 5n}]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Hermite polnomials
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[f,g,w,a,b,alpha,beta];
w[x_]:= Exp[-x^2];
a = -Infinity;
b = Infinity;

f[0,x_]=1;

alpha[0]= Integrate[x*w[x],{x,a,b}]/
          Integrate[ w[x],{x,a,b}];
f[1,x_]=x-alpha[0];

For[i=2,i<=5,i++,
    alpha[i-1]= 
      Integrate[Expand[x*f[i-1,x]*f[i-1,x]]*w[x],{x,a,b}]/
        Integrate[ Expand[f[i-1,x]*f[i-1,x]]*w[x],{x,a,b}];
    beta[i-1]= 
      Integrate[Expand[x*f[i-1,x]*f[i-2,x]]*w[x],{x,a,b}]/
        Integrate[Expand[ f[i-2,x]*f[i-2,x]]*w[x],{x,a,b}];
    f[i,x_]=Expand[(x-alpha[i-1])f[i-1,x]-beta[i-1] f[i-2,x]];
    ];
;[s]
3:0,0;4,1;19,0;565,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect]
MatrixForm[Table[Expand[f[n,x]/f[n,1]],{n,0,5}]]
MatrixForm[Table[Expand[HermiteH[n,x]/
HermiteH[n,1]],{n,0,5}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Nearly equioscillatory behavior
:[font = input; wordwrap; preserveAspect; endGroup]
n = 10;
Plot[{
  Sqrt[w[x]] HermiteH[n, x]}, 
  {x, -2Sqrt[n], 2Sqrt[n]},
  PlotRange -> 2(2n/E)^(n/2){-1, 1}]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Identity relating Hermite to Laguerre polynomials
:[font = input; wordwrap; preserveAspect; endGroup; endGroup; endGroup]
n = 5;

Print["H(", 2n, ",x) = "HermiteH[2n, x]];

constant = (-1)^(n)2^(2n) Factorial[n];
Expand[HermiteH[2n, x] - constant*
 LaguerreL[n, -1/2, x^2]]

Print["H(", 2n + 1, ",x) = "HermiteH[2n + 1, x]];

constant = (-1)^(n)2^(2n + 1) Factorial[n];
Expand[
HermiteH[2n + 1, x] - constant*
 x LaguerreL[n, 1/2, x^2]]
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Least square approximation
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Preliminaries
:[font = input; preserveAspect]
Off[General::spell1]
;[s]
3:0,0;4,1;19,0;21,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; preserveAspect]
Clear[LegendreApprox];
LegendreApprox[F_,n_] := Block[
{i,coef},
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           LegendreP[i,x]*F[x]],
        {x,-1,1}]/
        NIntegrate[Expand[
           LegendreP[i,x]*LegendreP[i,x]],
        {x,-1,1}];
    ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = input; preserveAspect]
Clear[ChebyshevUApprox];
ChebyshevUApprox[F_,n_] := Block[
{i,coef},
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           ChebyshevU[i,x]*F[x]*Sqrt[1-x^2]],
        {x,-1,1}]/
        NIntegrate[Expand[
           ChebyshevU[i,x]*ChebyshevU[i,x]*Sqrt[1-x^2]],
        {x,-1,1}];
    ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = input; preserveAspect]
Clear[ChebyshevTApprox];
ChebyshevTApprox[F_,n_] := Block[
{i,coef},
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           ChebyshevT[i,x]*F[x]/Sqrt[1-x^2]],
        {x,-1,1}]/
        NIntegrate[Expand[
           ChebyshevT[i,x]*ChebyshevT[i,x]/Sqrt[1-x^2]],
        {x,-1,1}];
    ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = input; preserveAspect]
Clear[LaguerreApprox];
LaguerreApprox[F_,n_,q_] := Block[
{i,coef},
If[q>0,
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           LaguerreL[i,q,x]*F[x]*x^q*Exp[-x]],
        {x,0,Infinity}]/
        NIntegrate[Expand[
           LaguerreL[i,q,x]*LaguerreL[i,q,x]*x^q*Exp[-x]],
        {x,0,Infinity}];
    ];,
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           LaguerreL[i,x]*F[x]*Exp[-x]],
        {x,0,Infinity}]/
        NIntegrate[Expand[
           LaguerreL[i,x]*LaguerreL[i,x]*Exp[-x]],
        {x,0,Infinity}];
    ];
  ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = input; preserveAspect; endGroup]
Clear[HermiteApprox];
HermiteApprox[F_,n_] := Block[
{i,coef},
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           HermiteH[i,x]*F[x]],
        {x,-Infinity,Infinity}]/
        NIntegrate[Expand[
           HermiteH[i,x]*HermiteH[i,x]],
        {x,-Infinity,Infinity}];
    ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Using Legendre polnomials of the FIRST kind
:[font = input; preserveAspect; startGroup]
Clear[F,a,b,p];
F[x_] = x Exp[-x]/Sqrt[1.1-x^2];
a = -1;
b = 1;
coef = LegendreApprox[F,5];
p[x_] := Expand[
  Sum[coef[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coef]-1}]]
p[x]
;[s]
3:0,0;16,1;49,0;176,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = input; preserveAspect; endGroup; endGroup]
Off[Plot::plnr]
Plot[{F[x],p[x]},{x,a,b}]
Plot[{F[x]-p[x]},{x,a,b}]
;[s]
3:0,0;4,1;14,0;68,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Approximation on FINITE  intervals [a,b] other than [-1,1]
:[font = text; inactive; preserveAspect; startGroup]
...using Legendre polnomials of the FIRST kind
:[font = input; preserveAspect; endGroup; endGroup]
Clear[F,a,b,p];
F[x_] = x Exp[-x];
a = -2;
b = 5;
f[x_]:= F[a+(x+1)/2(b-a)]


coef = LegendreApprox[f,5];
p[x_] := Expand[
  Sum[coef[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coef]-1}]]
P[x_] := p[2(x-a)/(b-a)-1]

Plot[{F[x],P[x]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]}]
Plot[{f[x],p[x]},{x,-1,1},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]}]
;[s]
3:0,0;16,1;35,0;362,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Notice that F(x)=f(x)*w(x),for nonnegative weight, then best to approxiamte f(x) by the 
p(x) using polynomials orthogonal with respect to the weight w(x).  The approximate in this case is w(x)*p(x)

This is illustrated below ...
;[s]
3:0,0;11,1;26,0;230,-1;
2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
...using Chebyshev polnomials of the SECOND kind
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
            Following is an empirical version of the function norm:
:[font = input; wordwrap; preserveAspect; endGroup]
Clear[funcNorm];
funcNorm[f_, p_, a_, b_] := If[p == Infinity,
    Max[Table[Abs[f[x]], {x, a, b, (b - a).001}]],
    NIntegrate[Abs[f[x]]^p, {x, a, b}]^(1/p)]
:[font = input; preserveAspect]
Off[Plot::plnr];
Clear[F,a,b,p,q,r];
F[x_] = x Exp[-x]Sqrt[1-x^2];
f[x_] = x Exp[-x];
a = -1;
b = 1;
n = 4;
;[s]
5:0,0;4,1;14,0;37,2;86,0;108,-1;
3:3,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use Legendre approximation of F[x], not taking weight into account
;[s]
7:0,1;4,2;12,1;30,0;34,1;36,3;39,1;67,-1;
4:1,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;1,13,10,Courier,1,12,65535,0,0;
:[font = input; preserveAspect; endGroup]
coefL = LegendreApprox[F,n];
p[x_] := Expand[
  Sum[coefL[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coefL]-1}]]
Print["Legendre approximation without any fuss..."]
Print["Thus p[x] is THE approximate of F[x]"]
Print["p[x]=",p[x]]
Clear[e];
Print[" "];
e[x_] := F[x]-p[x];
Print["The Infinity-norm of the error is ",NaiveInfNorm=funcNorm[e,Infinity,a,b]]
Print["The 2-norm of the error is ",NaiveTwoNorm=funcNorm[e,2,a,b]]

Plot[{F[x],p[x],
F[x]-q[x]Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed Legendre approximate"]

Plot[{F[x]-p[x]},{x,a,b},
PlotLabel->"Error for Legendre fit"]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use ChebyshevU approximation of f[x], optimally taking weight into account
;[s]
7:0,1;4,2;14,1;32,0;36,1;38,0;47,1;75,-1;
3:2,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;
:[font = input; preserveAspect; endGroup]
coefU = ChebyshevUApprox[f,n];
q[x_] := Expand[
  Sum[coefU[[i1+1]] ChebyshevU[i1,x],
  {i1,0,Length[coefU]-1}]]
Print["ChebyshevU approximation after dividing out weight..."]
Print["Thus w[x]q[x] is THE approximate of F[x]"]
Print["q[x]=",q[x]]
Print[" "];
Clear[e];
e[x_] := F[x]-q[x]Sqrt[1-x^2];
Print["The Infinity-norm of the error is ",OPTInfNorm=funcNorm[e,Infinity,a,b]]
Print["The 2-norm of the error is ",OPTTwoNorm=funcNorm[e,2,a,b]]

Plot[{F[x],q[x]Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed weighted ChebyshevU approximate"]

Plot[{F[x]-q[x]Sqrt[1-x^2]},{x,a,b},
PlotLabel->"Error for weighted ChebyshevU fit"]

;[s]
3:0,0;286,1;297,0;671,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use Legendre approximation of f[x], nonoptimally taking weight into account
;[s]
7:0,1;4,2;12,1;30,0;34,1;39,0;48,1;76,-1;
3:2,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;
:[font = input; preserveAspect; endGroup]
coefL = LegendreApprox[f,n];
r[x_] := Expand[
  Sum[coefL[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coefL]-1}]]
Print["Legendre approximation after dividing out weight..."]
Print["Thus w[x]r[x] is THE approximate of F[x]"]
Print["r[x]=",r[x]]
Print[" "];
Clear[e];
e[x_] := F[x]-r[x]Sqrt[1-x^2];
Print["The Infinity-norm of the error is ",NOPTInfNorm=funcNorm[e,Infinity,a,b]]
Print["The 2-norm of the error is ",NOPTTwoNorm=funcNorm[e,2,a,b]]

Plot[{F[x],r[x]Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed weighted Legendre approximate"]

Plot[{F[x]-r[x]Sqrt[1-x^2]},{x,a,b},
PlotLabel->"Error for weighted Legendre fit"]


;[s]
3:0,0;281,1;292,0;665,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Comparing the three approaches:
:[font = input; preserveAspect; endGroup; endGroup]
Print["Naive infinity-norm: ",NaiveInfNorm];
Print["Naive two-norm: ",NaiveTwoNorm];
Print["Optimal infinity-norm: ",OPTInfNorm];
Print["Optimal two-norm: ",OPTTwoNorm];
Print["Nonoptimal infinity-norm: ",NOPTInfNorm];
Print["Nonoptimal two-norm: ",NOPTTwoNorm];
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
...using Chebyshev polnomials of the FIRST kind
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
            Following is an empirical version of the function norm:
:[font = input; wordwrap; preserveAspect; endGroup]
Clear[funcNorm];
funcNorm[f_, p_, a_, b_] := If[p == Infinity,
    Max[Table[Abs[f[x]], {x, a, b, (b - a).001}]],
    NIntegrate[Abs[f[x]]^p, {x, a, b}]^(1/p)]
:[font = input; preserveAspect]
Off[Plot::plnr,Power::infy,NIntegrate::slwcon];
Clear[F,a,b,p,q,r];
F[x_] = x Exp[-x]/Sqrt[1-x^2];
f[x_] = x Exp[-x];
a = -1;
b = 1;
n = 4;
;[s]
5:0,0;4,1;45,0;68,2;118,0;140,-1;
3:3,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use Legendre approximation of F[x], not taking weight into account
;[s]
7:0,1;4,2;12,1;30,0;34,1;36,3;39,1;67,-1;
4:1,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;1,13,10,Courier,1,12,65535,0,0;
:[font = input; preserveAspect; endGroup]
coefL = LegendreApprox[F,n];
p[x_] := Expand[
  Sum[coefL[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coefL]-1}]]
Print["Legendre approximation without any fuss..."]
Print["Thus p[x] is THE approximate of F[x]"]
Print["p[x]=",p[x]]
Clear[e];
Print[" "];
e[x_] := F[x]-p[x];
Print["The Infinity-norm of the error is ",NaiveInfNorm=funcNorm[e,Infinity,a,b]]
Print["The 2-norm of the error is ",NaiveTwoNorm=funcNorm[e,2,a,b]]

Plot[{F[x],p[x],
F[x]-q[x]Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed Legendre approximate"]

Plot[{F[x]-p[x]},{x,a,b},
PlotLabel->"Error for Legendre fit"]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use ChebyshevT approximation of f[x], optimally taking weight into account
;[s]
7:0,1;4,2;14,1;32,0;36,1;38,0;47,1;75,-1;
3:2,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;
:[font = input; preserveAspect; endGroup]
coefT = ChebyshevTApprox[f,n];
q[x_] := Expand[
  Sum[coefT[[i1+1]] ChebyshevT[i1,x],
  {i1,0,Length[coefT]-1}]]
Print["ChebyshevU approximation after dividing out weight..."]
Print["Thus w[x]q[x] is THE approximate of F[x]"]
Print["q[x]=",q[x]]
Print[" "];
Clear[e];
e[x_] := (f[x]-q[x])/Sqrt[1-x^2];
Print["The Infinity-norm of the error is ",OPTInfNorm=funcNorm[e,Infinity,a+.00001,b-.00001]]
Print["The 2-norm of the error is ",OPTTwoNorm=funcNorm[e,2,a+.00001,b-.00001]]

Plot[{F[x],q[x]/Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed weighted ChebyshevT approximate"]

Plot[{F[x]-q[x]/Sqrt[1-x^2]},{x,a,b},
PlotLabel->"Error for weighted ChebyshevT fit"]

;[s]
3:0,0;289,1;300,0;704,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Use Legendre approximation of f[x], nonoptimally taking weight into account
;[s]
7:0,1;4,2;12,1;30,0;34,1;39,0;48,1;76,-1;
3:2,13,9,Times,0,12,0,0,0;4,13,9,Times,0,12,65535,0,0;1,13,10,Courier,1,12,0,0,0;
:[font = input; preserveAspect; endGroup]
coefL = LegendreApprox[f,n];
r[x_] := Expand[
  Sum[coefL[[i1+1]] LegendreP[i1,x],
  {i1,0,Length[coefL]-1}]]
Print["Legendre approximation after dividing out weight..."]
Print["Thus w[x]r[x] is THE approximate of F[x]"]
Print["r[x]=",r[x]]
Print[" "];
Clear[e];
e[x_] := (f[x]-r[x])/Sqrt[1-x^2];
Print["The Infinity-norm of the error is ",NOPTInfNorm=funcNorm[e,Infinity,a+.00001,b-.00001]]
Print["The 2-norm of the error is ",NOPTTwoNorm=funcNorm[e,2,a+.00001,b-.00001]]

Plot[{F[x],r[x]/Sqrt[1-x^2]},{x,a,b},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Dashed weighted Legendre approximate"]

Plot[{F[x]-r[x]/Sqrt[1-x^2]},{x,a,b},
PlotLabel->"Error for weighted Legendre fit"]

;[s]
3:0,0;284,1;295,0;697,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Comparing the three approaches:
:[font = input; preserveAspect; endGroup; endGroup; endGroup]
Print["Naive infinity-norm: ",NaiveInfNorm];
Print["Naive two-norm: ",NaiveTwoNorm];
Print["Optimal infinity-norm: ",OPTInfNorm];
Print["Optimal two-norm: ",OPTTwoNorm];
Print["Nonoptimal infinity-norm: ",NOPTInfNorm];
Print["Nonoptimal two-norm: ",NOPTTwoNorm];
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Approximation on SEMI-FINITE intervals [a,infinity] 
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
...using Laguerre polnomials of the FIRST kind
:[font = input; preserveAspect; endGroup]
Clear[LaguerreApprox];
LaguerreApprox[F_,n_,q_] := Block[
{i,coef},
If[q>0,
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           LaguerreL[i,q,x]*F[x]*x^q*Exp[-x]],
        {x,0,200}]/
        NIntegrate[Expand[
           LaguerreL[i,q,x]*LaguerreL[i,q,x]*x^q*Exp[-x]],
        {x,0,200}];
    ];,
For[i=0,i<=n,i++,
coef[i]=
        NIntegrate[Expand[
           LaguerreL[i,x]*F[x]*Exp[-x]],
        {x,0,200}]/
        NIntegrate[Expand[
           LaguerreL[i,x]*LaguerreL[i,x]*Exp[-x]],
        {x,0,200}];
    ];
  ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Suppose F[x]~f[x]*Exp[-x]
;[s]
2:0,0;8,1;26,-1;
2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;
:[font = input; preserveAspect; endGroup; endGroup; endGroup]
Clear[F,a,b,p];
Off[General::under,NIntegrate::ncvb];
F[x_] = (x+2) Sin[x/2] Exp[-x];
a = -2;
b = Infinity;
q = 0;
n = 8;

Plot[{F[x]},{x,a,2 n},
PlotRange->1.8{-1,1},
PlotLabel->"Original F[x]"]


(* Translate F to [0,Infinity], and divide by weight to get f *)
f[x_]:= F[x+a]*Exp[x+a]
Plot[{f[x]},{x,a,4 n},
PlotRange->4n{-1,1},
PlotLabel->"Translated and unweighted f[x]"]

(* Approximate new f by p using Laguerre polynomials *)
coef = LaguerreApprox[f,n,q];
If[q==0,
  p[x_] := Expand[
    Sum[coef[[i1+1]] LaguerreL[i1,x],
    {i1,0,Length[coef]-1}]],
  p[x_] := Expand[
    Sum[coef[[i1+1]] LaguerreL[i1,q,x],
    {i1,0,Length[coef]-1}]]
  ];
Print["p[x]=",p[x]]

(* Translate p to define polynomial P *)
P[x_] := p[x-a]
Print["P[x]=",P[x]]

Plot[{f[x],p[x]},{x,0,4n},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotRange->4n{-1,1},
PlotLabel->"Dashed p[x] approximates f[x]"]

Plot[{f[x]-p[x]},{x,0,4n},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"f[x]-p[x]"]

Plot[{F[x],P[x]Exp[-x]},{x,a,2n},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotRange->1.8{-1,1},
PlotLabel->"Dashed Exp[-x]*P[x] approximates F[x]"]

Plot[{F[x]-P[x]Exp[-x]},{x,a,2n},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
PlotLabel->"Exp[-x]*P[x]-F[x]"]
;[s]
13:0,0;20,2;51,0;54,1;86,0;198,1;263,0;286,1;287,0;377,1;433,0;671,1;712,0;1261,-1;
3:7,7,10,Courier,1,12,0,0,0;5,7,10,Courier,1,12,65535,0,0;1,7,10,Courier,0,12,0,0,0;
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Least square approximation on finite set
;[s]
3:0,0;30,1;36,0;41,-1;
2:2,19,14,Times,1,18,0,0,0;1,19,14,Times,1,18,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Uniformly distributed points
:[font = text; inactive; preserveAspect; startGroup]
					Best fit on finite set of m+1 points
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[g,w,a,b,alpha,beta];
w[x_]:= 1;

m=10;
X = Table[-1+2(i1-1)/m,{i1,m+1}];


g[0,x_]=1;

alpha[0]= Sum[X[[i1]]*w[X[[i1]]],{i1,m+1}]/
          Sum[w[X[[i1]]],{i1,m+1}];
g[1,x_]=x-alpha[0];
 
For[i=2,i<=10,i++,
  alpha[i-1]=Sum[X[[i1]]*g[i-1,X[[i1]]]*g[i-1,X[[i1]]]*w[X[[i1]]],{i1,m+1}]/Sum[g[i-1,X[[i1]]]*g[i-1,X[[i1]]]*w[X[[i1]]],{i1,m+1}];
  
beta[i-1]=Sum[X[[i1]]*g[i-1,X[[i1]]]*g[i-2,X[[i1]]]*w[X[[i1]]],{i1,m+1}]/Sum[g[i-2,X[[i1]]]*g[i-2,X[[i1]]]*w[X[[i1]]],{i1,m+1}];

g[i,x_]=Expand[(x-alpha[i-1])g[i-1,x]-beta[i-1]g[i-2,x]];
];
;[s]
3:0,0;4,1;19,0;562,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect; endGroup]
MatrixForm[Table[Expand[g[n,x]],{n,0,5}]]
:[font = text; inactive; preserveAspect; startGroup]
					Example in book, page 56
:[font = input; preserveAspect]
Clear[LegendreFinite];
LegendreFinite[F_,n_] := Block[
{i,coef},
For[i=0,i<=n,i++,
coef[i]=
        Sum[
           g[i,X[[i1]]]*F[X[[i1]]],
        {i1,m+1}]/
        Sum[
           g[i,X[[i1]]]*g[i,X[[i1]]],
        {i1,m+1}];
    ];
Return[Table[coef[i1],{i1,0,n}]];
]
:[font = input; preserveAspect; endGroup; endGroup]
Clear[F,p,q];
F[x_] := Abs[x];
coef = LegendreFinite[F,8];
p[x_]:= Expand[Sum[coef[[i1]]*g[i1-1,x],
  {i1,Length[coef]}]];
p[t]

Plot[{F[x],p[x]},{x,-1,1},
PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]},
Prolog->{PointSize[.02],
Table[Point[{X[[i1]],p[X[[i1]]]}],{i1,Length[X]}]}]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Arcsine distributed points
:[font = text; inactive; preserveAspect; startGroup]
					Best fit on finite set of m+1 points
:[font = input; wordwrap; preserveAspect]
Off[General::spell1];
Clear[g,w,a,b,alpha,beta];
w[x_]:= 1;

m=10;
pi = N[Pi];
X = Table[Cos[(i1-1)pi/m],{i1,m+1}];


g[0,x_]=1;

alpha[0]= Sum[X[[i1]]*w[X[[i1]]],{i1,m+1}]/
          Sum[w[X[[i1]]],{i1,m+1}];
g[1,x_]=x-alpha[0];
 
For[i=2,i<=10,i++,
  alpha[i-1]=Sum[X[[i1]]*g[i-1,X[[i1]]]*g[i-1,X[[i1]]]*w[X[[i1]]],{i1,m+1}]/Sum[g[i-1,X[[i1]]]*g[i-1,X[[i1]]]*w[X[[i1]]],{i1,m+1}];
  
beta[i-1]=Sum[X[[i1]]*g[i-1,X[[i1]]]*g[i-2,X[[i1]]]*w[X[[i1]]],{i1,m+1}]/Sum[g[i-2,X[[i1]]]*g[i-2,X[[i1]]]*w[X[[i1]]],{i1,m+1}];

g[i,x_]=Expand[(x-alpha[i-1])g[i-1,x]-beta[i-1]g[i-2,x]];
];
;[s]
3:0,0;4,1;19,0;577,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;
:[font = input; wordwrap; preserveAspect; endGroup; endGroup]
MatrixForm[Table[Expand[g[n,x]],{n,0,5}]]
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
Usual least squares on discrete points
:[font = text; inactive; preserveAspect]
   If think of the error in fitting m points as   E(f)=sum [f(xi)-p(xi)]^2,
   then E(f) is a function with n+1 unknowns.  To find the minimum, we consider 
   the system Ax=b, where the rows of Ax are p(xi), and the rows of b are f(xi).  Multiplying 
   both sides by A', the transpose of A, yields the normal equations (A'A)x=A'b. the solution
   to this system is the minimum of the least squares E(f) above.
;[s]
3:0,0;304,1;310,0;412,-1;
2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;
:[font = input; preserveAspect]
m=10;
n=8;
Clear[p,q];
:[font = text; inactive; Cclosed; preserveAspect; startGroup]
           Approximate on uniform nodes
:[font = input; preserveAspect; endGroup]
X = Table[-1+2(i1-1)/m,{i1,m+1}];

mat = Transpose[Table[Table[X[[i1]]^j1,{j1,n}],{i1,m+1}]];
mat = Transpose[Prepend[mat,Table[1,{m+1}]]];
dat = Table[{F[X[[i1]]]},{i1,m+1}];
ata = Transpose[mat].mat;
atb = Transpose[mat].dat;
Ucoef = Flatten[Inverse[ata].atb];

q[x_] := Sum[Ucoef[[i1]] x^(i1-1),{i1,n+1}]

:[font = text; inactive; Cclosed; preserveAspect; startGroup]
           Approximate on Chebyshev nodes
:[font = input; preserveAspect; endGroup]
pi = N[Pi];
X = Table[Cos[(i1-1)pi/m],{i1,m+1}];

mat = Transpose[Table[Table[X[[i1]]^j1,{j1,n}],{i1,m+1}]];
mat = Transpose[Prepend[mat,Table[1,{m+1}]]];
dat = Table[{F[X[[i1]]]},{i1,m+1}];
ata = Transpose[mat].mat;
atb = Transpose[mat].dat;
CCoef = Flatten[Inverse[ata].atb];

p[x_] := Sum[CCoef[[i1]] x^(i1-1),{i1,n+1}]

:[font = input; preserveAspect; endGroup; endGroup]
p[t]
q[t]
Plot[{F[t],q[t],p[t]},{t,-1,1},
PlotStyle->{Dashing[{1,0}],Dashing[{.05,.05}],Dashing[{.01,.02}]}]
Plot[{F[t]-q[t],F[t]-p[t]},{t,-1,1},
PlotStyle->{Dashing[{.05,.05}],Dashing[{.01,.02}]}]
^*)