/************************************************** * Evaluate a Legendre Polynomial P(x) of Order n * * for argument x by using Horner's Rule * * ----------------------------------------------- * * SAMPLE RUN: * * * * give order: 7 * * * * give argument x: 0.5 * * * * P(x) = 0.223145 * * * * ----------------------------------------------- * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * **************************************************/ #include #include #define NMAX 12 #define TRUE 1 /***************************************************** * Legendre series coefficients evaluation subroutine * * by means of recursion relation. The order of the * * polynomial is n. The coefficients are returned in * * A(i) by increasing powers. * *****************************************************/ void Legendre_Coeff(int n, double *A) { int i,j; double B[NMAX][NMAX]; // Establish p0 and p1 coefficients B[0][0]=1.0; B[1][0]=0.0; B[1][1]=1.0; // Return if order is less then 2 if (n > 1) { for (i=2; i<=n; i++) { B[i][0] = -(i-1)*B[i-2][0]/i; for (j=1; j<=i; j++) { // Basic recursion relation B[i][j]=(i+i-1)*B[i-1][j-1]-(i-1)*B[i-2][j]; B[i][j] /= i; } } for (i=0; i<=n; i++) A[i]=B[n][i]; } } void HORNER (int N, double *A, double X0, int K, bool B, double *R) { /*-------------------------------------------------------------------- ! EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD ! ! INPUTS: ! N ORDER OF POLYNOMIAL ! A VECTOR OF SIZE N+1 STORING THE COEFFICIENTS OF ! POLYNOMIAL IN DECREASING ORDERS OF POWERS ! I.E. P(X) = A(1)*X**N+A(2)*X**(N-1)+...+A(N)*X+A(N+1) ! X0 GIVEN ARGUMENT ! K MAXIMUM ORDER OF DERIVATIVES ASKED FOR ! B FLAG ! = .TRUE. EVALUATION ONLY ! = .FALSE. DETERMINATION OF POLYNOMIAL COEFFICIENTS ! NEAR X0 ! R VECTOR OF SIZE K+1 CONTAINS: ! - VALUES P(X0), P'(X0), P''(X0),..PK(X0), ! IF B IS TRUE AND IF K < N ! - THE N+1 COEFFICIENTS OF POLYNOMIAL ! Q(X-X0) = R(1)+R(2)*(X-X0)+...+R(N+1)*(X-X0)**(N), ! IF B IS FALSE AND IF K = N. ! ! REFERENCE: ! ALGORITHM 337, COLLECTED ALGORITHMS FROM CACM, W.PANKIEWICKZ !-------------------------------------------------------------------*/ double RR; int I,J,L,NMJ; RR = A[1]; for (I=1; I<=K+1; I++) R[I] = RR; for (J=2; J<=N+1; J++) { R[1] = R[1]*X0 + A[J]; NMJ = N-J+1; if (NMJ > K) L = K; else L = NMJ; for (I=2; I<=L+1; I++) R[I] = R[I]*X0 + R[I-1]; } if (B) { L = 1; for (I=2; I<=K+1; I++) { L *= (I-1); R[I] *= L; } } } // Evaluate P(x) of order n for argument x double Eval(int n, double *A, double x) { int i; double A1[NMAX], R[NMAX]; for (i=n; i>=0; i--) A1[n-i+1] = A[i]; HORNER(n,A1,x,0,TRUE,R); return R[1]; } void main() { double A[NMAX]; int n; double x; printf("\n give order: "); scanf("%d", &n); printf("\n give argument x: "); scanf("%lf", &x); Legendre_Coeff(n, A); printf("\n P(x) = %f\n\n", Eval(n,A,x)); } // end of file eval_leg.cpp