Changeset 3514 in Sophya for trunk/SophyaLib/NTools

Timestamp:

Aug 8, 2008, 4:17:50 PM (17 years ago)

Author:

cmv

Message:

3.14159 _> M_PI + commentaires dans nbgauleg cmv 08/08/2008

File:

• trunk/SophyaLib/NTools/nbmath.c (modified) (2 diffs)

trunk/SophyaLib/NTools/nbmath.c

@@ -1121 +1121 @@
 |  w[] = valeur des coefficients associes
 | REMARQUES:
-|  - x et w doivent au moins etre dimensionner a n.
-|  - l'integration est rigoureuse si sur l'intervalle d'integration
-|    la fonction f(x) peut etre approximee par un polynome
-|    de degre 2*m (monome le + haut x**(2*n-1) )
+|  - x et w doivent au moins etre dimensionnes a n.
+|  - l'integration est exacte sur l'intervalle [x1,x2]
+|    si f(x) est un polynome de degre inferieur ou egal a 2*n-1
 |  - Voir la fonction Integ_Fun pour un calcul d'ordre 8
 --
@@ -1131 +1130 @@
 {
    int m,j,i;
-   double z1,z,xm,xl,pp,p3,p2,p1;
+   double z1,z,xm,xl,pp,p3,p2,p1; // Need high precision
 
    m=(n+1)/2;
-   xm=0.5*(x2+x1);
-   xl=0.5*(x2-x1);
-   for (i=1;i<=m;i++)  {
-      z=cos(3.141592654*(i-0.25)/(n+0.5));
-      do {
+   xm=0.5*(x2+x1);     // The root are symetric in the intervalle
+   xl=0.5*(x2-x1);     // we just need to find half of them
+   for (i=1;i<=m;i++)  {   // Loop over the desired root
+     z=cos(M_PI*(i-0.25)/(n+0.5));  // Starting with the approx for the i-ieme root
+     do {       // Entering the main loop of refinement by Newton's method
          p1=1.0;
          p2=0.0;
-         for (j=1;j<=n;j++) {
+         for (j=1;j<=n;j++) { // Loop up the reccurence relation to get the Legendre polynomial at z
             p3=p2;
             p2=p1;
             p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
          }
+         // p1 is now the desired Legendre polynomial. We next compute pp, its derivative,
+         // by a standard relation involving also p2, the polynomial of one lower order
          pp=n*(z*p1-p2)/(z*z-1.0);
          z1=z;
-         z=z1-p1/pp;
+         z=z1-p1/pp;   // Newton's method
       } while (fabs(z-z1) > EPS_gauleg);
-      x[i-1]=xm-xl*z;
-      x[n-i]=xm+xl*z;
-      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
-      w[n-i]=w[i-1];
+      x[i-1]=xm-xl*z;                   // Scale the root to the desired interval,
+      x[n-i]=xm+xl*z;                   // and put in its symetric counterpart.
+      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);  // Compute the weight
+      w[n-i]=w[i-1];                    // and its symetric counterpart
    }
 }