4
$\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$
  • $\begingroup$ Thanks David and @georhe2079 I'm quite new here :-) $\endgroup$ – AstroFloor Nov 9 '16 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 Nov 9 '16 at 21:45
  • $\begingroup$ just tried it, NIntegrate[ .. WorkingPrecision -> 100, MaxRecursion -> 50] gives all real>0 roots. $\endgroup$ – george2079 Nov 9 '16 at 21:49
  • $\begingroup$ @george2079, brilliant! That definitely solved the problem! Thank you so much! $\endgroup$ – AstroFloor Nov 9 '16 at 22:38
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$
  • $\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 Nov 19 '16 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 Nov 19 '16 at 15:04
  • $\begingroup$ Sorry, I don't know how to ask you this in a better way Thanks in advance $\endgroup$ – AstroFloor Nov 19 '16 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 Nov 19 '16 at 15:38
  • $\begingroup$ @AstroFloor I added some explanation. HTH $\endgroup$ – Michael E2 Nov 19 '16 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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