Thread Subject:

What is the maxumum number of abscissas could be used for the Gauss-Hermite quadrature?

When I integrate an infinite quadrature from -inf to inf. I use Gauss-hermite Quadrature rule. First of all, I need to calculate the abscissas and weights. But I found I could only obtain up to 388 point of weights and abscissas using Jacobi matrix as follows,
function [x,w] = Ngqzero(n)
d = sqrt((1:(n-1))/2);
T = diag(d,1)+diag(d,-1);
[V,D] = eig(T);
w = V(1,:).^2;
w = w*sqrt(pi);
x = diag(D);
[x,ind] = sort(x);
w = w(ind);
w = w';
end

For which when n<=388, the weights are symmetric and continuously increase from 0 to the two ends. However, when n>388, the weights at the two end of the weights vector become 0.Which is unreasonable obviously because the weights should become larger at those points instead of 0.
So could anyone know some other program to obtained correct abscissas and weights more than 388, because I need large abscissas range to calculate the quadrature.
Thank you !

"Ala " <yingzhangning@126.com> wrote in message <hmkfdu$gu4$1@fred.mathworks.com>...
> When I integrate an infinite quadrature from -inf to inf. I use Gauss-hermite Quadrature rule. First of all, I need to calculate the abscissas and weights. But I found I could only obtain up to 388 point of weights and abscissas using Jacobi matrix as follows,
> function [x,w] = Ngqzero(n)
> d = sqrt((1:(n-1))/2);
> T = diag(d,1)+diag(d,-1);
> [V,D] = eig(T);
> w = V(1,:).^2;
> w = w*sqrt(pi);
> x = diag(D);
> [x,ind] = sort(x);
> w = w(ind);
> w = w';
> end
>
> For which when n<=388, the weights are symmetric and continuously increase from 0 to the two ends. However, when n>388, the weights at the two end of the weights vector become 0.Which is unreasonable obviously because the weights should become larger at those points instead of 0.
> So could anyone know some other program to obtained correct abscissas and weights more than 388, because I need large abscissas range to calculate the quadrature.
> Thank you !
----------------
  Ala, I ran your program for n = 7 and obtained this result:

[x,w] =
  -2.65196135683523 0.00097178124510
  -1.67355162876747 0.05451558281913
  -0.81628788285896 0.42560725261013
  -0.00000000000000 0.81026461755681
   0.81628788285896 0.42560725261013
   1.67355162876747 0.05451558281913
   2.65196135683523 0.00097178124510

  The weights in the second column are those that would approximate the integral

 int(-inf,inf), exp(-x^2)*f(x) dx, by the sum

 w(1)*f(x1)+w(2)*f(x2) + ... + w(7)*f(x7)

according to the Gauss-Hermite quadrature rule. As you see, this kind of weight does not increase towards either end but rather decreases and becomes very small towards those ends.

  To obtain an approximation to an integral of the form

 int(-inf,inf) f(x) dx

each of the w values should be multiplied by the corresponding exp(x(i)^2). In that case the weight values would indeed increase from the middle (though not from zero) towards either end, as you have described.

  This suggests the possible source of your trouble for n > 388. I am guessing that for such large values of n, your original w values have underflowed to an exact zero value at the two ends if you obtained them using the matlab code you showed. In spite of a subsequent multiplication by exp(x(i)^2) to obtain the second kind of weight, such a value would stubbornly remain at a precise zero.

  Remember that the very smallest possible positive nonzero quantity that an IEEE double precision number can possess is 2^(-1074). There is nothing between that and an exact zero. (Moreover the numbers between 2^(-1022) and 2^(-1074), which are called denormalized numbers, gradually lose precision as their size decreases and will yield increasingly inaccurate results.)

  I would suggest that your problem should be reformulated so as to make use of the first of the two above concepts of weighting rather than the second one. However you do it, for large values of n the contributions from those end terms in your approximating sum should be quite small or else the Gauss-Hermite approximation will not accurate.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hms4bl$nb1$1@fred.mathworks.com>...
> .......
> This suggests the possible source of your trouble for n > 388. I am guessing that for such large values of n, your original w values have underflowed to an exact zero value at the two ends if you obtained them using the matlab code you showed.
> ......

  Friday in this thread I guessed that for very large n your difficulties might be due to a hypothesized underflow to zero in the endpoint w values obtained from the eigenvectors in your computation. However upon reflection, I think an underflow would be highly unlikely in this situation. What I think actually happens is that the theoretical values of these endpoint w's would be so small as to become completely obscured by the natural rounding errors inherent in the eigenvector computation, and the values you obtain at these ends would then be totally meaningless. Furthermore, I suspect this happens long before you get to n = 388. At n = 32 the end values of w are down to around 7.3e-23 and in the process they have already lost six decimal places of accuracy out of the usual sixteen.

  In my opinion, to use values of n in the hundreds you need to have far higher accuracy than could possibly be provided by 64-bit double precision numbers no matter what algorithm is used. I would suggest that you try making appropriate use of the Matlab symbolic toolbox with a very much larger number of digits of precision setting for computing the quadrature abscissas and weights.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hms4bl$nb1$1@fred.mathworks.com>...
> "Ala " <yingzhangning@126.com> wrote in message <hmkfdu$gu4$1@fred.mathworks.com>...
> > When I integrate an infinite quadrature from -inf to inf. I use Gauss-hermite Quadrature rule. First of all, I need to calculate the abscissas and weights. But I found I could only obtain up to 388 point of weights and abscissas using Jacobi matrix as follows,
> > function [x,w] = Ngqzero(n)
> > d = sqrt((1:(n-1))/2);
> > T = diag(d,1)+diag(d,-1);
> > [V,D] = eig(T);
> > w = V(1,:).^2;
> > w = w*sqrt(pi);
> > x = diag(D);
> > [x,ind] = sort(x);
> > w = w(ind);
> > w = w';
> > end
> >
> > For which when n<=388, the weights are symmetric and continuously increase from 0 to the two ends. However, when n>388, the weights at the two end of the weights vector become 0.Which is unreasonable obviously because the weights should become larger at those points instead of 0.
> > So could anyone know some other program to obtained correct abscissas and weights more than 388, because I need large abscissas range to calculate the quadrature.
> > Thank you !
> ----------------
> Ala, I ran your program for n = 7 and obtained this result:
>
> [x,w] =
> -2.65196135683523 0.00097178124510
> -1.67355162876747 0.05451558281913
> -0.81628788285896 0.42560725261013
> -0.00000000000000 0.81026461755681
> 0.81628788285896 0.42560725261013
> 1.67355162876747 0.05451558281913
> 2.65196135683523 0.00097178124510
>
> The weights in the second column are those that would approximate the integral
>
> int(-inf,inf), exp(-x^2)*f(x) dx, by the sum
>
> w(1)*f(x1)+w(2)*f(x2) + ... + w(7)*f(x7)
>
> according to the Gauss-Hermite quadrature rule. As you see, this kind of weight does not increase towards either end but rather decreases and becomes very small towards those ends.
>
> To obtain an approximation to an integral of the form
>
> int(-inf,inf) f(x) dx
>
> each of the w values should be multiplied by the corresponding exp(x(i)^2). In that case the weight values would indeed increase from the middle (though not from zero) towards either end, as you have described.
>
> This suggests the possible source of your trouble for n > 388. I am guessing that for such large values of n, your original w values have underflowed to an exact zero value at the two ends if you obtained them using the matlab code you showed. In spite of a subsequent multiplication by exp(x(i)^2) to obtain the second kind of weight, such a value would stubbornly remain at a precise zero.
>
> Remember that the very smallest possible positive nonzero quantity that an IEEE double precision number can possess is 2^(-1074). There is nothing between that and an exact zero. (Moreover the numbers between 2^(-1022) and 2^(-1074), which are called denormalized numbers, gradually lose precision as their size decreases and will yield increasingly inaccurate results.)
>
> I would suggest that your problem should be reformulated so as to make use of the first of the two above concepts of weighting rather than the second one. However you do it, for large values of n the contributions from those end terms in your approximating sum should be quite small or else the Gauss-Hermite approximation will not accurate.
>
> Roger Stafford



Yes, the smallest weight (w(1)) for N = 388 is 1e-323 about 2^(-1073) and if N >388, as you predicted, w(1)<1e-323 and overflow occurs. But the problem is, even such a large N cannot make the integrals as low as enough to do the truncation if we ignore the precision.
Another thing I need to mention is,
Then I try to decrease the error limit through changing the significant digits to make it 32 digits (which is the largest digits MATLAB can use) instead of 16 digits which is the default one.
I found when N = 42, w0 (1) = 1e-31 (the first element of the original weights), which means even for the 32 digits, the largest N which is precise for use is 42. While when I compare the weights w1 (0) to that obtained using double precision mode (16digit), they are equal for the first 16 digits, which seems mean that the default 16digits MATLAB using is not such a bad one (but who knows how bad?), the same thing happens when N =112, the first 16 digits are the same even if both modes are not so accurate.
I am thinking of a way to produce the second weight directly without resorting to obtaining the first weights first, in order to produce the right weights when N>388.
Because for the second weight, it is always relative large.

"Ala " <yingzhangning@126.com> wrote in message <hua1ln$a6t$1@fred.mathworks.com>...

> Then I try to decrease the error limit through changing the significant digits to make it 32 digits (which is the largest digits MATLAB can use) instead of 16 digits which is the default one.

No, you cannot change the precision in matlab to
32 digits. This is not an option, unless you are using
the symbolic toolbox, in which case you can set any
precision you desire.

So what version of matlab are you using? Perhaps a
release from the future?

John

Yes, I use the symbolic box to define the precision (32digits).
But as I mentioned, even for 32digits, when N = 42, w0(1) = 1e-31, which means the weights and abscissas when N >42may be inaccurate.
So under this circumstance, how could I using the Gauss-Hermit integration rule for an infinite integration which needs N is over 388?
Could anyone help me with that? Thank you very much!

Btw, I am using Matlab7.7.0(R2008b).

"Ala " <yingzhangning@126.com> wrote in message <hvuq5k$jiu$1@fred.mathworks.com>...
> Yes, I use the symbolic box to define the precision (32digits).
> But as I mentioned, even for 32digits, when N = 42, w0(1) = 1e-31, which means the weights and abscissas when N >42may be inaccurate.
> So under this circumstance, how could I using the Gauss-Hermit integration rule for an infinite integration which needs N is over 388?
> Could anyone help me with that? Thank you very much!
>
> Btw, I am using Matlab7.7.0(R2008b).
- - - - - - - - -
  What is preventing you from setting the symbolic toolbox precision to much higher values than 32 decimal digits? Why not 100, 400, 1000, whatever you need to achieve the accuracy you are after? I have often set the symbolic precision on my symbolic toolbox to well over 100 and it does just fine. Where did you get the idea that 32 was a limit?

Roger Stafford