(***********************************************************************)
(*                                                                     *)
(* Roman Maurer <roman.maurer@fmf.uni-lj.si>, 30-Jul-1996              *)
(* 1. domača naloga iz Numerične analize:                              *)
(* Naloga 7. c: Aproksimacija funkcije sin(5x) na intervalu [0,pi/2]   *)
(* s polinomom stopnje največ 10:                                      *)
(*   - Bernsteinov polinom                                             *)
(*   - polinom najboljše enakomerne aproksimacije                      *)
(*   - p.n.e.a po metodi najmanjših kvadratov                          *)
(*                                                                     *)
(***********************************************************************)

(* interval na katerem aproksimiramo funkcijo *)
a=0;
b=N[Pi/2]

(* stopnja aproksimacijskega polinoma *)
n = 10

(* r := residual *)
r[x_] := N[ f[x] - p[x] ];

(* - računamo (bisekcijo ipd.) na natančnost eps
   - delilne točke, oddaljene <eps2 imamo za enake *)
eps = 10^(-10);
eps2 = 10^(-3)

(* funkcija, ki jo aproksimiramo *)
f[x_] := Sin[5*x]


Sgn::usage = "Sgn[f,a,b] vrne True, če sta f[a] in f[b] istega predznaka."

Sgn[f_,a_,b_] := Sign[N[f[a]]] == Sign[N[f[b]]]


Bisection::usage = "Bisection[f,a,b] z bisekcijo poišče ničlo funkcije f
na intervalu [a,b]. Predpostavi, da je funkcija v krajiščih različno
predznačena."

Bisection[f_,left_,right_] := Module[ {c, a=left, b=right},
   While[ Abs[b-a] > eps,
      c = N[ (a+b)/2 ];
      If[ Sgn[f,c,a], a=c, b=c ]
   ];
   Return[c] ]


SecondNorm::usage = "SecondNorm[f] vrne drugo normo funkcije f."

SecondNorm[f_] := Sqrt[ NIntegrate[ (f[x])^2, {x,a,b} ] ]


FindMax::usage = "FindMax[f,a,b] poišče maksimalno vrednost funkcije |f|
na intervalu [a,b]. Predpostavi, da je f' v krajiščih nasprotno
predznačena in da je vmes en sam ekstrem."

FindMax[f_,a_,b_] := Abs[ f[ Bisection[f',a,b] ] ]


(***************************************************************************)

Bernstein::usage = "Bernstein[f,x] vrne n-ti Bernsteinov polinom
za funkcijo f v točki x."

Bernstein[f_,x_] :=  Simplify[  (1-x)^n * f[a] +
   Sum[ Binomial[n,i] * (1-x)^(n-i) * x^i * f[a+i*(b-a)/n], {i, 1, n} ]
   ] /. x->(x-a)/(b-a)

(***************************************************************************)


BestApproximationOnSet::usage = "BestApproximationOnSet[f,list] izračuna
polinom najboljše enakomerne aproksimacije na množici delilnih točk
iz seznama list. Vrne seznam {m,a_0,a_1,...,a_(n+1)}, v katerem je
prvi element m napaka aproksimacije.
Ostali elementi seznama določajo koeficiente p.n.e.a."

BestApproximationOnSet[f_,pts_] := N[ LinearSolve[
   (* matrika = (1,-1,1,-1,1,...) + Van der Mondova matrika *)
   Transpose[ Join[ {Table[ (-1)^i, {i, 0, n+1} ]},
      {Table[ 1, {n+2} ]},
      Table[ pts[[i]]^j, {j, 1, n}, {i, 1, n+2} ] ] ],
   (* desna stran je vektor vrednosti funkcije f v točkah *)
   N[ f /@ pts ] ] ]


MaxArg::usage = "MaxArg[list] vrne zaporedno številko največjega
elementa v seznamu list."

MaxArg[l_List] := 0 /; l == {}

MaxArg[l_List] := Module[ {i, maxl, indmaxl, len=Length[l]},
   maxl = l[[1]];  (* največja vrednost *)
   indmaxl=1;      (* indeks maksimalnega elementa *)
   i=2;            (* indeks trenutnega elementa *)
   While[ i<=len,
      If[ l[[i]] > maxl,
         maxl = l[[i]];
         indmaxl = i
      ];
      i++
   ];
   Return[indmaxl]
] /; l =!= {}


Exchange::usage = "Exchange[r,points,y] opravi Remesovo zamenjavo točke y
s točko iz seznama points glede na residual r."

Exchange[r_,pt_List,y_] := Module[ {i, newp},
   (* i := indeks delilne točke, levo od y, ali 0, če je y < pt[[1]] *)
   i = Count[ Negative[pt-y], True ];
   If[ (i>=1) && (i<=n+1), (*then*)
      If[ Sgn[r,pt[[i]],y], (*then*)
         newp = Union[ {y}, Delete[pt,i] ],
      (*else*) newp = Union[ {y}, Delete[pt,i+1] ]
      ],
   (*else*) If[ i==0, (*then*)
      If[ Sgn[r,pt[[1]],y], (*then*)
         newp = Union[ {y}, Delete[pt,1] ],
      (*else*) newp = Union[ {y}, Delete[pt,n+2] ]
      ],
   (*else*) If[ Sgn[r,pt[[n+2]],y], (*then*)
      newp = Union[ {y}, Delete[pt,n+2] ],
   (*else*) newp = Union[ {y}, Delete[pt,1] ]
   ] ] ];
   Return[newp]
]


RemesOneStep::usage = "RemesOneStep[f,list] naredi en korak Remesovega prvega
postopka za funkcijo f na točkah iz seznama list. Vrne novi seznam točk.
Spremeni globalno funkcijo p, ki zdaj določa aproksimacijski polinom."

RemesOneStep[f_,pts_] := Module[ {rZero, rExtreme, coeffs, d, y, newpts, p, r},
   (* poiščemo p.n.a. na delilnih točkah *)
   coeffs = BestApproximationOnSet[f,pts];
   (* p := aproksimacijski polinom *)
   p[x_] = coeffs[[2]] + Sum[ coeffs[[i]] * x^(i-2), {i, 3, n+2} ];
   r[x_] = N[ f[x]-p[x] ];
   rZero = Table[ Bisection[ r, pts[[i]], pts[[i+1]] ], {i, 1, n+1} ];
   rExtreme = Join[ {a},
      Table[ Bisection[ r', rZero[[i]], rZero[[i+1]] ], {i, 1, n} ],
      {b} ];
   (* y := globalni maksimum funkcije r(x) na [a,b] *)
   y = rExtreme[[ MaxArg[ Abs[r[rExtreme]] ] ]];
   (* d := razlika residuala v točki y in v delilnih točkah *)
   d = Abs[ r[y] ] - Abs[ coeffs[[1]] ];
   If [ (d>eps), (*then*)
      newpts = Exchange[r,pts,y],
   (*else*) newpts = pts ];
   If [ Max[Abs[newpts-pts]] < eps2, (*then*)
      Return[pts],    (* točke oddaljene <eps2 so za nas enake *)
   (*else*) Return[newpts] ]
]


Remes::usage = "Remes[f,x] vrne polinom najboljše enakomerne aproksimacije
za funkcijo f v točki x, dobljen po prvem Remesovem postopku."

Remes[f_,x_] := Module[
{points = N[ Table[ a + (b-a)*j/(n+1), {j, 0, (n+1)} ] ], coeffs},
(*   Remove[p]; *)
   (* points := začetne delilne točke *)
   (* opravimo korak Remesovega postopka, dokler gre *)
   (* polinom je skrit v seznamu koeficientov *)
   coeffs = BestApproximationOnSet[f,FixedPoint[RemesOneStep[f,#]&,points]];
   Return[ coeffs[[2]] + Sum[ coeffs[[i]] * x^(i-2), {i, 3, n+2} ] ]
]

(***************************************************************************)

ScalarProduct::usage = "ScalarProduct[f,g] vrne diskretni skalarni produkt
funkcij f in g."

ScalarProduct[f_,g_] := Sum[ f[x]*g[x] /. x->a+(b-a)*j/(n+1), {j, 0, (n+1)} ]


DiscreteMinSqr::usage = "DiscreteMinSqr[f,x] vrne polinom najboljše
aproksimacije po diskretni metodi najmanjših kvadratov."

DiscreteMinSqr[f_,x_] := Module[ {g, G},
   (* g[i] := i-ti bazni vektor *)
   g[1] = 1&;
   g[j_] := #^(j-1)& /; (j>1) && (j<=n+1);
   G = N[Table[ ScalarProduct[g[i],g[j]], {i,1,n+1}, {j,1,n+1} ]];
   coeffs = LinearSolve[G, Table[ ScalarProduct[f,g[i]], {i,1,n+1}]];
   Return[ coeffs[[1]] + Sum[ coeffs[[i]] * x^(i-1), {i,2,n+1}] ] ]