From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: What is the maxumum number of abscissas could be used for the Gauss-Hermite quadrature?
Date: Fri, 4 Jun 2010 05:07:03 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 60
Message-ID: <hua1ln$a6t$>
References: <hmkfdu$gu4$> <hms4bl$nb1$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1275628023 10461 (4 Jun 2010 05:07:03 GMT)
NNTP-Posting-Date: Fri, 4 Jun 2010 05:07:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1904048
Xref: comp.soft-sys.matlab:642093

"Roger Stafford" <> wrote in message <hms4bl$nb1$>...
> "Ala " <> wrote in message <hmkfdu$gu4$>...
> > 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.