5
$\begingroup$

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]$ :-)

$\endgroup$
4
  • $\begingroup$ Thanks David and @georhe2079 I'm quite new here :-) $\endgroup$
    – AstroFloor
    Commented Nov 9, 2016 at 21:28
  • $\begingroup$ 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 $\endgroup$
    – george2079
    Commented Nov 9, 2016 at 21:45
  • $\begingroup$ just tried it, NIntegrate[ .. WorkingPrecision -> 100, MaxRecursion -> 50] gives all real>0 roots. $\endgroup$
    – george2079
    Commented Nov 9, 2016 at 21:49
  • $\begingroup$ @george2079, brilliant! That definitely solved the problem! Thank you so much! $\endgroup$
    – AstroFloor
    Commented Nov 9, 2016 at 22:38

1 Answer 1

3
$\begingroup$

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
Quadraturepoints =
   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
Quadraturepoints2 =
   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_].

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – AstroFloor
    Commented Nov 19, 2016 at 15:00
  • $\begingroup$ 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? $\endgroup$
    – AstroFloor
    Commented Nov 19, 2016 at 15:04
  • $\begingroup$ Sorry, I don't know how to ask you this in a better way Thanks in advance $\endgroup$
    – AstroFloor
    Commented Nov 19, 2016 at 15:13
  • $\begingroup$ 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.) $\endgroup$
    – Michael E2
    Commented Nov 19, 2016 at 15:38
  • $\begingroup$ @AstroFloor I added some explanation. HTH $\endgroup$
    – Michael E2
    Commented Nov 19, 2016 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.