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