# Why do I get quadrature points $\notin [0,\infty]$, when making own Gaussian quadrature rule with range $[0,\infty]$?

For my project, I'm trying to calculate in Mathematica my own Gaussian quadrature rule (i.e. the quadrature points and weights) for a certain model with a parameter $x$ with the following (Maxwallian) probability density function: $\rho(x) = \sqrt{\frac{2}{\pi}}\frac{1}{265^3} x^2 e^{-x^2 / (2\times 265^2)}$ on range $[0,\infty]$

Let ${P}(N)$ be the space with all polynomials of degree less or equal to $N$ and let $\{ \psi_n \}_{n=o}^{n=N}$ form an orthogonal basis with respect to our $\rho(x)$ on $[0,\infty]$. The $N$ Gaussian quadrature points are then defined as the roots of the $N$-th basis function, $\psi_N$. (Thus, to get the quadrature points one only has to determine the $N$-th orthogonal basis function and find the $N$ distinct roots.) This strategy seems to work for my other probability density functions in my project.

However, when I try to solve the quadrature points for this $\rho(x)$ I have the problem that when deriving more than $N>7$ quadrature points, the NSolve gives me negative and sometimes non-identical quadrature points as a solution, and I don't know why. I've tried a few things and think that or the Orthogonalization I do is wrong for basis functions of higher degree than 6, or maybe the way I use NSolve does not work in this case, or maybe something else? I would be happy with any suggestions to solve this.

Here is my Mathematica code:

Orthogonalbasis =  Orthogonalize[x^Range[0, 10],
NIntegrate[#1*#2* Sqrt[(2/\[Pi])]*(1/(265^3)) *  x^2 * Exp[(-x^2/(2(265)^2))],
{x, 0, Infinity}, MaxRecursion -> 16] &];
Quadraturepoints = x /. NSolve[Orthogonalbasis[[8]], x ]


{-5723.14, 106.843, 264.616, 462.599, 692.796, 957.759, 1279.71}

(Where I have used the polynomial basis $\{1,x,...,x^{10} \}$, and used MaxRecursion to solve another problem).

But I don't understand why I get -5723.14 as quadrature point, which is not in $[0,\infty]$. Even imaginary points occur when trying to solve the roots of higher degree basis functions, for instance solving the 8th and 9th degree polynomial basis function gives:

x /. NSolve[Orthogonalbasis[[9]], x]


{-1439.52 - 770.095 I, -1439.52 + 770.095 I, 109.61, 269.795, 469.023, 699.113, 962.56, 1281.01}

x /. NSolve[Orthogonalbasis[[10]], x]


{-1496.94 - 1528.1 I, -1496.94 + 1528.1 I, 107.32, 265.383, 463.198, 692.617, 956.015, 1275.05, 4776.29}

Any suggestions on how I can derive $N$ quadrature points in the range $[0,\infty]$? Thank you!

NB: these are actually the half range generalized Hermite polynomials and points on $[0,\infty]$ :-)

• Thanks David and @georhe2079 I'm quite new here :-) – AstroFloor Nov 9 '16 at 21:28
• I think its a precision problem. Have you tried Integrate? It should work in theory but make take a long time. You might also work with higher working precision in NIntegrate – george2079 Nov 9 '16 at 21:45
• just tried it, NIntegrate[ .. WorkingPrecision -> 100, MaxRecursion -> 50] gives all real>0 roots. – george2079 Nov 9 '16 at 21:49
• @george2079, brilliant! That definitely solved the problem! Thank you so much! – AstroFloor Nov 9 '16 at 22:38

You can use Integrate to give a general, symbolic formula for the inner product and use N[] to speed up the orthogonalization. One needs to use fairly high precision because the coefficients in the symbolic result become quite large in the OP's example.

(* OP's approach, with george2079's WorkingPrecision tip, for comparison *)
Orthogonalbasis = Orthogonalize[x^Range[0, 10],
NIntegrate[#1*#2*Sqrt[(2/π)]*(1/(265^3))*x^2*Exp[(-x^2/(2 (265)^2))],
{x, 0, Infinity}, MaxRecursion -> 16,
PrecisionGoal -> 8, WorkingPrecision -> 32] &]; // AbsoluteTiming
x /. NSolve[Orthogonalbasis[[8]], x, WorkingPrecision -> MachinePrecision]
(*
{2.77581, Null}
{87.3421, 219.779, 389.161, 587.07, 811.184, 1066.68, 1375.93}
*)

(* inner product *)
ClearAll[ip];
SetAttributes[ip, Listable];
ip[n_Integer, s_] =             (* basic integral, with parameter s *)
Integrate[x^n*Sqrt[(2/π)]*(1/(s^3))*Exp[(-x^2/(2 (s)^2))],
{x, 0, Infinity}, Assumptions -> n > 0 && n ∈ Integers && s > 0];
ip[f_, g_, s_] :=               (* inner product of f, g, with parameter s *)
#.ip[Range[2, Length@# + 1], s] &@ CoefficientList[f*g, x]

(* orthogonalization, with s = 265, as in the OP *)
Orthogonalbasis = Orthogonalize[x^Range[0, 10], N[ip[##, 265], 100] &]; // AbsoluteTiming
x /. NSolve[Orthogonalbasis[[8]], x, WorkingPrecision -> MachinePrecision]
(*
{0.046954, Null}
{87.3421, 219.779, 389.161, 587.07, 811.184, 1066.68, 1375.93}
*)


The ip[] method is about 65 times faster than the NIntegrate[] approach to orthogonalize the degree-10 basis.

Addendum: How ip computes the inner product

There are two components, ip[n, s], which represents the integral of the monomial x^n times the function Sqrt[(2/π)]*(1/(s^3))*Exp[(-x^2/(2 (s)^2))], and ip[f, g, s], which represents the inner product of polynomials f and g (assumed to depend on x). The inner product has a weight function Sqrt[(2/π)]x^2(1/(s^3))*Exp[(-x^2/(2 (s)^2))], which isx^2times the "reduced" weight function inip[n, s].

I defined ip to have the attribute Listable so that ip[{2, 3, 4, 5}, s] would return a list of the values of the integrals of {x^2, x^3, x^4, x^5}. These can be used to evaluate the integral of a polynomial by taking the dot product of these values with the coefficient list of the polynomial.

For instance for the inner product of a linear and quadratic polynomial, we would integrate

f*g   -> a0 + a1 x + a2 x^2 + a3 x^3


with respect to the inner product weight function, or equivalently the polynomial f * g * x^2 with respect to the reduced weight function in ip[n, s]. The degrees of the monomials are {2, 3, 4, 5}. The coefficient list is given by

CoefficientList[f*g, x]   -> {a0, a1, a2, a3}


The list of degrees is given by

Range[2, Length@CoefficientList[f*g, x] + 1]   -> {2, 3, 4, 5}


Thus the inner product is given by

CoefficientList[f*g, x] . ip[Range[2, Length@CoefficientList[f*g, x] + 1], s]


which is coded more succinctly (and more efficiently) in the definition of ip[f_, g_, s_].

• Thank you for your helpful answer. However, I think I do not fully understand your code: When I try to use your method on finding the orthogonal system with respect to the weight function \rho(x) = x^(-2.3) on range [1,300], I'm doing something wrong. – AstroFloor Nov 19 '16 at 15:00
• I'm using your code, but now for x^(-2.3) instead of the maxwallian: ip[n_Integer, s_] =(I don't need to use s, so have set s=1 *) Integrate[x^n*(s^3 * x^(-2.3)), {x, 1, 300}, Assumptions -> n > 0 && n [Element] Integers && s > 0]; ip[f_, g_, s_] :=#.ip[Range[2, Length@# + 1], s] &@CoefficientList[fg, x] Orthogonalbasistest = Orthogonalize[x^Range[0, 11], N[ip[##, 1], 1000] &] But when I test it, I don;t get 0, e.g.: Integrate[ Orthogonalbasistest[[1]]*Orthogonalbasistest[[3]]* x^(-2.3), {x, 1, 300}] (* 0.015378 *) Has this something to do with the s? – AstroFloor Nov 19 '16 at 15:04
• Sorry, I don't know how to ask you this in a better way Thanks in advance – AstroFloor Nov 19 '16 at 15:13
• Why are you multiplying the product of the basis elements by x^2.3 instead of the weight function in the inner product? And why are you integrating only from 1 to 300, instead of 0 to Infinity`? (I don't see any relationship between your test integral and the inner product in your question. I'm not sure why the test integral ought to be zero.) – Michael E2 Nov 19 '16 at 15:38
• @AstroFloor I added some explanation. HTH – Michael E2 Nov 19 '16 at 15:57