(* opolynom.wl -- Gram-Schmidt to Find Orthogonal Polynomials *)

dp=Function[{p,q},Integrate[p*q,{x,-1,1}]]
nm=Function[p,Sqrt[dp[p,p]]]

(* Standard polynomial basis *)

u1=1
u2=x
u3=x^2
u4=x^3
u5=x^4

(* Perform Gram-Schmidt algorithm *)

tv1=u1
v1=tv1/nm[tv1]
tv2=u2-dp[u2,v1]*v1
v2=tv2/nm[tv2]
tv3=u3-dp[u3,v1]*v1-dp[u3,v2]*v2
v3=tv3/nm[tv3]
tv4=u4-dp[u4,v1]*v1-dp[u4,v2]*v2-dp[u4,v3]*v3
v4=tv4/nm[tv4]
tv5=u5-dp[u5,v1]*v1-dp[u5,v2]*v2-
      dp[u5,v3]*v3-dp[u5,v4]*v4
v5=tv5/nm[tv5]
Simplify[v5]

(* Find the roots of the 4th degree polynomial *)

S=Solve[v5==0,x]
x1=N[x/.S[[1]],15]
x2=N[x/.S[[2]],15]
x3=N[x/.S[[3]],15]
x4=N[x/.S[[4]],15]

(* Also compute an exact integral for comparison *)

exact=Integrate[Log[x+2],{x,-1,1}]
N[exact,15]
