Mathematics algorithm

[dude_1967]

Hello,

I know that there are some forum members with considerable mathematics knowledge. So I'll just ask...

Can anyone please explain exactly how to compute weights and abscissae on the integral range -1 to +1 for N-degree Gaussian-Quadrature?

I have read the relevant chapter in Numerical Recipes and I studied the Wolfram Research page, but for some reason I just didn't get it. The numbers always just appear in tabulated form and for some reason I never understand how they got there.

I really want to be able to understand the exact numerical formulae for these numbers without resorting to using tables or Mathematica.

Thanks for any help.

Sincerely, Chris.

[Philip Nicoletti]

maybe this will help (around page 20 in the PDF file ...)

http://www.nek.uu.se/Utbildning/Fosk...g/MEnotes1.pdf

[galathaea]

If you are needing these types of results, you are most likely doing some kind of numerical analysis, so first let me suggest the book Introduction to Numerical Analysis by F. B. Hildebrand. This is not to suggest that you need an introductory book or even that the book is really good. I merely suggest it because its Dover price (at least when I bought it) was a suggested US $14.95 for 669 pages, and I have found that it is usually a good source for those simple little explanations that I couldn't find elsewhere, so I think it is a good value (particularly since you can probably find it used for 4 to 8 bucks US pretty easily).

That said, its section on Gaussian quadrature is really good. Now, you do not mention a weight function, so I will assume you are discussing a true Gaussian and not a Christoffel, weighted form. For the interval [-1, 1], things get very easy. Summarising what is needed for your calculation, the answer is

Integral([-1,1], f(x) dx) = Sum(1 <= k <= m, Hk f(xk)) + error

where

Hk = 2 / (m P[m - 1; xk] P'[m; xk])

P[n; x] is the nth Legendre polynomial

and xk is the kth root of P[m; x].

There is an example given for m = 3, where the 3rd Legendre polynomial is

P[3; x] = (1 / 2) (5 x^3 - 3 x)

whose roots are therefore

x1 = - &radic;(15) / 5
x2 = 0
x3 = &radic;(15) / 5

using also the 2nd Legebdre

P[2; x] = (1 / 2) (3 x^2 - 1)

and taking the root of the third (the P' in the formula), we can calculate the Hk to be

H1 = 5/9
H2 = 8/9
H3 = 5/9

So basically, the calculation is first of the Legendre polynomials, the roots of the Nth for N-degree quadrature giving the abscissas, and the Hk (a. k. a. Christoffel numbers) given by a straightforward formula on other Legendre polynomials (and / or a derivative -- the above formulation being one of many different equivalent forms due to all the transformation rules the polynomials obey) at the abscissas.

Just for completeness, the above simple calculation gives the final form for the quadrature:

Integral([-1, 1], f(x) dx) =
(1/9) [5 f(-&radic;(15)/5) + 8 f(0) + 5 f (&radic;(15)/5)] + error

[dude_1967]

Right on!

The link from Philip looks relevant. In addition, the worked-out example from galathaea is just what I was looking for.

You see, I am just putting the wraps on a new extended precision floating point library. I've been working on this kind of stuff for a few years now. I've always tested similar libraries using calculations of simple elementary functions and some special functions such as bessel or gamma functions.

(BTW: Yes, I have Mathematica and Matlab, but I just love doing the extended precision numerical stuff and learning alongside these mighty programs).

This time, I want to gain some experience with extended precision numerical integration. When I was studying in college, I used a lot of stuff like Simpson's-type methods and they seemed OK. However, for the exptended precision stuff, I think that it really makes sense to cut down on the number of function evaluations---thus the interest in Gaussian-Quadrature. However, I just didn't get it!

Now I can go on.

Thanks a bunch. This will keep me ticking for a few weeks.

Sincerely, Chris.

P.S. I was asking about true Guass-Quad in the range -1 to +1 with a weighting function of 1, so the example is completely relevant.