(***********************************************************************)
(* NA2.M                                                               *)
(* Roman Maurer <roman.maurer@fmf.uni-lj.si>, 11-Aug-1996              *)
(* 2. domača naloga iz Numerične analize:                              *)
(* Naloga 5. f: Numerično določi uteži in vozle Gaussovih              *)
(* integracijskih pravil:                                              *)
(* $\int_0^1 f(t) \sqrt(x) dx \approx \sum_{i=1}^n w_i f(x_i)$         *)
(* za n=2,3,...,10                                                     *)
(*                                                                     *)
(***********************************************************************)


(* interval na katerem računamo integral *)
l=0;
r=1

(* utež *)
W::usage = "W[x] je utež prostora funkcij."

W[x_] = Sqrt[x]


(* n=2,3,...,10 *)
maxn=10


eps = 10^(-10);

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] ]


PrFrm::usage = "PrFrm[x] daje prvih 8 decimalk števila x kot niz."

PrFrm[x_] := StringTake[ToString[N[Round[10^8*x]/10^8,8]] <> "00000000", 10]


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

Print["Roman Maurer, uporabna matematika, 11-Aug-1996"];
Print["2. domača naloga iz Numerične analize:"];
Print["Naloga 5. f: Numerično določi uteži in vozle Gaussovih"];
Print["integracijskih pravil na intervalu [0,1]:"];
Print["integral f(t) sqrt(x) dx =cca. w_1 f(t_1) + ... + w_n f(t_n)"];
Print["za n=2,3,...,10"];


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

ScalarProduct[f_,g_] := Integrate[ W[x]*f[x]*g[x], {x,l,r} ]


SqrNorm::usage = "SqrNorm[n] vrne kvadrat norme n-tega ortogonalnega
polinoma p[n,x]."

SqrNorm[-1] = 0;
SqrNorm[0] = Integrate[ W[x], {x,l,r}]

p::usage = "p[i,x] vrne ortogonalni (glede na utež W) polinom stopnje i."

b[0]=0;
a[0]= ScalarProduct[ #&, 1& ] / SqrNorm[0];

p[-1,x_] = 0;
p[0,x_] = 1;
p[1,x_] = x-a[0];

knot::usage = "knot[i] vrne seznam vozlov Gaussovega integracijskega
pravila z i členi."



knot[1] = {Bisection[p[1,#]&,l,r]}


For[n=2, n<=maxn, n++,
   Print[ " " ];
   Print["\tn=",n];
   SqrNorm[n-1] = ScalarProduct[p[n-1,#]&,p[n-1,#]&];
   a[n-1] = ScalarProduct[ #*p[n-1,#]&, p[n-1,#]& ] / SqrNorm[n-1];
   b[n-1] = SqrNorm[n-1] / SqrNorm[n-2];
   p[n,x_] =  Simplify[ (x-a[n-1])*p[n-1,x] - b[n-1]*p[n-2,x] ];

   (* pri izračunu vozlov upoštevamo, da gre za ničle OG polinomov in
      da se ničle zaporednih polinomov prepletajo *)
   knot[n] = {Bisection[ p[n,#]&,l,knot[n-1][[1]] ]};
   For[i=1,i<n-1,i++,
      knot[n] = Union[ knot[n],
         {Bisection[ p[n,#]&,knot[n-1][[i]],knot[n-1][[i+1]] ]} ];
   ];
   knot[n] = Union[ knot[n],
      {Bisection[ p[n,#]&,knot[n-1][[n-1]],r ]} ];

   Print[ "   t_i \t           w_i" ];
   (* določimo uteži *)
   For[i=1,i<=n,i++,
      w[n,i] = SqrNorm[n-1] / (p[n-1,x]*D[p[n,x],x]) /. x->knot[n][[i]];
      Print[ PrFrm[ knot[n][[i]] ], " \t", PrFrm[ w[n,i] ] ]
   ]
]


GaussIntegrateWithSqrt::usage = "GaussIntegrateWithSqrt[f] izračuna
integral funkcije f(x)*sqrt(x) na intervalu (0,1) po Gaussovem
kvadraturnem pravilu na 10 točkah."

GaussIntegrateWithSqrt[ f_ ] :=
   Sum[ w[maxn,i] * f[ knot[maxn][[i]] ], {i,1,maxn}];


Compare::usage = "Compare[f[x],x] primerja Gaussov približek
(z desetimi točkami) za integral f(x)*sqrt(x) s pravo vrednostjo."

Compare[f_,x_] := Module[ {h, gauss, int},
   h[t_] = f /. x->t;
   gauss = GaussIntegrateWithSqrt[h];
   int = NIntegrate[ h[t] Sqrt[t], {t,0,1} ];
   Print["Gauss\tprava vrednost"];
   Print[ gauss, "\t", int];
   Print["napaka=", N[ gauss-int,3] ];
   Print["relativna napaka=", 100*N[gauss/int-1,5],"%" ];
]