(**********************************************
  Approximation algorithms 
 **********************************************)

(* Question 3  *)
(* phi : function, n : no. of basis functions
   x,f : sample points (size m) *)
LSQ[phi_,n_,x_,f_] := 
 Block[{m = Length[x],A,A1,i,j,b,c},
       A = Table[phi[x[[i]],j],{i,m},{j,n}]; (* A : mxn *)
       A1 = Transpose[A].A;
       b = Transpose[A].f;
       c = LinearSolve[A1,b];
       Return[c];];

LSError[F_,f_,x_] := Block[{error=0.0,i,n=Length[x]},
			   Do[error=error+(f[[i]]-F[x[[i]]])^2,{i,n}];
			   Return[Sqrt[error]];];




x = {-1.,-.5,-.25,0.,.25,.5,1.};
f = {4.1,2.9,2.1,1.4,.91,.3,-.7};

(* 3a *)
Phi3[x_,i_] := Return[Switch[i,1,1,2,x,3,Exp[x]]];
c3 = LSQ[Phi3,3,x,f];
F3[x_] := Return[c3[[1]]+c3[[2]]*x+c3[[3]]*Exp[x]];
P1 = Plot[F3[x],{x,-1.0,1.0}];
LSError[F3,f,x]
(* F3 = 1.05992 - 2.92505*x + 0.427231*e^x *)
(* Least Square Error = 0.1583 *)


(* 3b *)
Phi2[x_,i_] := Return[Switch[i,1,1,2,x]];
c2 = LSQ[Phi2,2,x,f];
F2[x_] := Return[c2[[1]]+c2[[2]]*x];
P2 = Plot[F2[x],{x,-1.,1.}];
LSError[F2,f,x]
(* F2 = 1.57286 - 2.43714*x *)
(* Least Square Error = 0.298867 *)

(* Question 2 *)

(* Define Inner Product *)
Pb2 = N[Pi/2.0];
GProduct[g_,h_] := Return[Integrate[g*h,{x,-Pb2,Pb2}]];

(* Use this to get the orthogonal polynomials *)
Pm1[x_] := 0.;
P0[x_] := 1.;
P1[x_] := x;
P2[x_] := -0.822467 + x^2;
P3[x_] := -1.48044*x + x^3;
P4[x_] := 0.521834 - 2.11491*x^2  + x^4;
P5[x_] := Pnext[P3[x],P4[x]];

(* Use the polynomials to get the coeffecients *)
Coef[Pn_] := GProduct[Sin[x],Pn[x]]/GProduct[Pn[x],Pn[x]];
Pnext[Pim_,Pi_] := Return[Simplify[(x-GProduct[x*Pi,Pi]/GProduct[Pi,Pi])*Pi
			  -GProduct[Pi,Pi]/GProduct[Pim,Pim]*Pim]];


(* The least square polynomial hence obtained is as follows *)
(* Coef[P0] = 0, Coef[P1] = 0.774037, Coef[P2] = 0.0;
   
  Answer : Poly = 0.774037*x; *)

(* Justification : We want to minimize, the integral of square of the error,
   The integral corresponds to a inner product Integral(P(x)*Sin(x)).
   The least value is obtained by mapping Sin(x) orthogonally onto
   the space of functions defined by the polynomial of order 2.

   The result is obtained by a linear combination of all orthogonal
   polynomials of upto degree 2 (over -Pi/2 to Pi/2), where the coeffecients
   of the polynomials are given by the orthogonal projection of Sin(x) onto
   the polynomial *)


   
