Changeset 20021

Timestamp:

01/29/16 17:15:13 (9 years ago)

Author:

Eric.Larour

Message:

CHG: replace legendre.cpp by a better routine.

Location:

issm/trunk-jpl/src/c/shared/Numerics

Files:

• legendre.cpp (modified) (1 diff)
• numerics.h (modified) (1 diff)

issm/trunk-jpl/src/c/shared/Numerics/legendre.cpp

@@ -137 +137 @@
 #include "./types.h"
 #include "../Exceptions/exceptions.h"
+#include "../MemOps/MemOps.h"
 #include "./recast.h"
 
-IssmDouble legendre(IssmDouble Pn1, IssmDouble Pn2, IssmDouble x, int n){
-
-        if (n<0)_error_("legendre error message: order required should be >0");
-
-        if(n==0)return 1;
-        if(n==1)return x;
-
-        return ( (2*n-1)*x*Pn1 - (n-1)*Pn2 ) /n;
+double *p_polynomial_value ( int m, int n, IssmDouble* x){
+
+        /******************************************************************************{{{/
+        Purpose:
+
+        P_POLYNOMIAL_VALUE evaluates the Legendre polynomials P(n,x).
+
+        Discussion:
+
+        P(n,1) = 1.
+        P(n,-1) = (-1)^N.
+        | P(n,x) | <= 1 in [-1,1].
+
+        The N zeroes of P(n,x) are the abscissas used for Gauss-Legendre
+        quadrature of the integral of a function F(X) with weight function 1
+        over the interval [-1,1].
+
+        The Legendre polynomials are orthogonal under the inner product defined
+        as integration from -1 to 1:
+
+        Integral ( -1 <= X <= 1 ) P(I,X) * P(J,X) dX
+        = 0 if I =/= J
+        = 2 / ( 2*I+1 ) if I = J.
+
+        Except for P(0,X), the integral of P(I,X) from -1 to 1 is 0.
+
+        A function F(X) defined on [-1,1] may be approximated by the series
+        C0*P(0,x) + C1*P(1,x) + ... + CN*P(n,x)
+        where
+        C(I) = (2*I+1)/(2) * Integral ( -1 <= X <= 1 ) F(X) P(I,x) dx.
+
+        The formula is:
+
+        P(n,x) = (1/2^N) * sum ( 0 <= M <= N/2 ) C(N,M) C(2N-2M,N) X^(N-2*M)
+
+        Differential equation:
+
+        (1-X*X) * P(n,x)'' - 2 * X * P(n,x)' + N * (N+1) = 0
+
+        First terms:
+
+        P( 0,x) =      1
+        P( 1,x) =      1 X
+        P( 2,x) = (    3 X^2 -       1)/2
+        P( 3,x) = (    5 X^3 -     3 X)/2
+        P( 4,x) = (   35 X^4 -    30 X^2 +     3)/8
+        P( 5,x) = (   63 X^5 -    70 X^3 +    15 X)/8
+        P( 6,x) = (  231 X^6 -   315 X^4 +   105 X^2 -     5)/16
+        P( 7,x) = (  429 X^7 -   693 X^5 +   315 X^3 -    35 X)/16
+        P( 8,x) = ( 6435 X^8 - 12012 X^6 +  6930 X^4 -  1260 X^2 +   35)/128
+        P( 9,x) = (12155 X^9 - 25740 X^7 + 18018 X^5 -  4620 X^3 +  315 X)/128
+        P(10,x) = (46189 X^10-109395 X^8 + 90090 X^6 - 30030 X^4 + 3465 X^2-63)/256
+
+        Recursion:
+
+        P(0,x) = 1
+        P(1,x) = x
+        P(n,x) = ( (2*n-1)*x*P(n-1,x)-(n-1)*P(n-2,x) ) / n
+
+        P'(0,x) = 0
+        P'(1,x) = 1
+        P'(N,x) = ( (2*N-1)*(P(N-1,x)+X*P'(N-1,x)-(N-1)*P'(N-2,x) ) / N
+
+        Licensing:
+
+        This code is distributed under the GNU LGPL license.
+
+        Modified:
+
+        08 August 2013
+
+        Author:
+
+                John Burkardt
+
+                        Reference:
+
+                        Milton Abramowitz, Irene Stegun,
+                                   Handbook of Mathematical Functions,
+                                   National Bureau of Standards, 1964,
+                                   ISBN: 0-486-61272-4,
+                                   LC: QA47.A34.
+
+                                           Daniel Zwillinger, editor,
+                                   CRC Standard Mathematical Tables and Formulae,
+                                   30th Edition,
+                                   CRC Press, 1996.
+
+                                           Parameters:
+
+                                           Input, int M, the number of evaluation points.
+
+                                           Input, int N, the highest order polynomial to evaluate.
+                                           Note that polynomials 0 through N will be evaluated.
+
+                                           Input, double X[M], the evaluation points.
+
+                                           Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre
+                                           polynomials of order 0 through N.
+        }}}*/
         
+        int i;
+        int j;
+        IssmDouble* v=NULL;
+
+        if ( n < 0 ) return NULL;
+
+        v = xNew<IssmDouble>(m*(n+1));
+
+        for ( i = 0; i < m; i++ ) v[i+0*m] = 1.0;
+        if ( n < 1 ) return v;
+
+        for ( i = 0; i < m; i++ ) v[i+1*m] = x[i];
+
+        for ( j = 2; j <= n; j++ ) {
+                for ( i = 0; i < m; i++ ) {
+                        v[i+j*m] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[i+(j-1)*m]
+                                        - ( IssmDouble ) (     j - 1 ) *        v[i+(j-2)*m] )
+                                / ( IssmDouble ) (     j     );
+                }
+        }
+
+        return v;
 }

issm/trunk-jpl/src/c/shared/Numerics/numerics.h

@@ -34 +34 @@
 void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
 int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d,IssmDouble X[3], int *num);
-IssmDouble  legendre(IssmDouble Pn1, IssmDouble Pn2, IssmDouble x, int i);
+IssmDouble  legendre(IssmDouble Pn1, IssmDouble Pn2, IssmDouble x, int n);
+IssmDouble*  p_polynomial_value ( int m, int n, IssmDouble* x);
 
 int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm);