Integrating with Gaussian quadratures method. Precision and Stability

776 Views Asked by At

We have gaussian quadratures method: $$ \int^\infty_0 f(x) dx \approx \sum^n_iw_i f(x_i) $$

I want to compute abscissas and the weights for various $n$ ( where $n$ will be users input) for lognormal distribution function. $$f(x) = \frac{1}{\sqrt{2\pi}x\sigma} e^{-\frac{(log(x) - \mu)^2}{2\sigma^2}}$$

Analytical integration over the $R$ gives $\int_R f = 1$. This is my sanity check that Gaussian quadratures works fine.

  1. Why does higher $n$ gives worse result ? ( Is it some kind of over-fitting?)
  2. Can I get the same precision for different $n$? ( I want my result be stable)

The best result I get is aprox. $0.99978$ (for $n \approx 20$ I get it from Laguerre quadrature method. )

I ask this because I saw black box realization of Gauss quadrature for lognormal function, where the sum of weigths multiplied by function of abscissas always give sum equal to 1 (for any $n$ I put inside).


I did ask question before: Gaussian quadratures problems with precision and stability

There I got a better understanding but still not clear picture what is wrong in what I am doing.

I also found that the error should be (for Legendre) $$\frac{(b-a)^{2n+1} (n!)^4}{(2n+1) (2n!)^3} f^{(2n)} (\xi) , \quad \xi \in [a,b]$$ so I believe that black box is doing smth. different, as no visible error in sum of abscissas and weights is present.

  1. Is it right? If yes, what shall I do to get similar effect as in the black box?

Also I have tried to split interval of integration into two parts: $$[0,a] -\text{integrate by Legendre, } [a,\infty) -\text{integrate by Laguerre, }$$ My result: $0.999992$ (21 points for each method,i.e. in total 41 distinct points) - which is obviously better than it was before but I still want to get maximum precision and stability possible. I read a Davis, Rabinowitz book "Methods of numerical integration", which was the most rigorous about the topic and yet I understood something :), but it did not cover the precision, stability as I wanted.

There were propositions to use Stroud,Secrest method, Konrad Patterson, preassign abscissas but I do not know still how to do that in my case.

P.S. I do need to have as a solution 2 vectors: abscissas and weights that give me the sum of $\sum_i w_i f(\text{abscissas}_i)=1 $.

When I have obtained the $0.999992$ value, I used Legendre method for interval $[0,a]$ and for $[a,\infty)$ interval I have used Laguerre method. Of course, for Laguerre I plugged in the integral extra exponent $$\int^\infty_a e^{-x}\underset{g(x)}{(e^x f(x))}dx \approx \sum^n_i w_i g(x_i)$$

I have been asked to do precisely Gaussian quadratures of lognormal as it gives $(x_i,w_i)$ that will be used for other analysis further. Anyway, I myself can not understand how that black box works ( I do not know how it works, but in its description it says: plug mean, standard deviation and get 2 vectors of abscissas and weights), and I doubt myself if it is doing the same stuff as I have been asked.

That is why I ask questions (1)-(3), I want to understand what are the limits in precision, stability and if am I doing everything correct. I am computing that approximate integral, as black box gives 1.

I am sorry if I am not too precise ( if I can improve smth -please tell)