Ice Sheet System Model  4.18
Code documentation
Functions
legendre.cpp File Reference

compute legendre polynomial at specific x. This is inspired from the p_polynomial_value package, written by John Burkardt (see below for full documentation). Here, we apply the function differently: we assume previous two evaluations are provided, and we compute the level n evaluation: More...

#include "./types.h"
#include "../Exceptions/exceptions.h"
#include "../MemOps/MemOps.h"
#include "./recast.h"

Go to the source code of this file.

Functions

IssmDoublep_polynomial_value (int m, int n, IssmDouble *x)
 

Detailed Description

compute legendre polynomial at specific x. This is inspired from the p_polynomial_value package, written by John Burkardt (see below for full documentation). Here, we apply the function differently: we assume previous two evaluations are provided, and we compute the level n evaluation:

Definition in file legendre.cpp.

Function Documentation

◆ p_polynomial_value()

IssmDouble* p_polynomial_value ( int  m,
int  n,
IssmDouble x 
)

Definition at line 142 of file legendre.cpp.

142  {
143 
144  /******************************************************************************{{{/
145  Purpose:
146 
147  P_POLYNOMIAL_VALUE evaluates the Legendre polynomials P(n,x).
148 
149  Discussion:
150 
151  P(n,1) = 1.
152  P(n,-1) = (-1)^N.
153  | P(n,x) | <= 1 in [-1,1].
154 
155  The N zeroes of P(n,x) are the abscissas used for Gauss-Legendre
156  quadrature of the integral of a function F(X) with weight function 1
157  over the interval [-1,1].
158 
159  The Legendre polynomials are orthogonal under the inner product defined
160  as integration from -1 to 1:
161 
162  Integral ( -1 <= X <= 1 ) P(I,X) * P(J,X) dX
163  = 0 if I =/= J
164  = 2 / ( 2*I+1 ) if I = J.
165 
166  Except for P(0,X), the integral of P(I,X) from -1 to 1 is 0.
167 
168  A function F(X) defined on [-1,1] may be approximated by the series
169  C0*P(0,x) + C1*P(1,x) + ... + CN*P(n,x)
170  where
171  C(I) = (2*I+1)/(2) * Integral ( -1 <= X <= 1 ) F(X) P(I,x) dx.
172 
173  The formula is:
174 
175  P(n,x) = (1/2^N) * sum ( 0 <= M <= N/2 ) C(N,M) C(2N-2M,N) X^(N-2*M)
176 
177  Differential equation:
178 
179  (1-X*X) * P(n,x)'' - 2 * X * P(n,x)' + N * (N+1) = 0
180 
181  First terms:
182 
183  P( 0,x) = 1
184  P( 1,x) = 1 X
185  P( 2,x) = ( 3 X^2 - 1)/2
186  P( 3,x) = ( 5 X^3 - 3 X)/2
187  P( 4,x) = ( 35 X^4 - 30 X^2 + 3)/8
188  P( 5,x) = ( 63 X^5 - 70 X^3 + 15 X)/8
189  P( 6,x) = ( 231 X^6 - 315 X^4 + 105 X^2 - 5)/16
190  P( 7,x) = ( 429 X^7 - 693 X^5 + 315 X^3 - 35 X)/16
191  P( 8,x) = ( 6435 X^8 - 12012 X^6 + 6930 X^4 - 1260 X^2 + 35)/128
192  P( 9,x) = (12155 X^9 - 25740 X^7 + 18018 X^5 - 4620 X^3 + 315 X)/128
193  P(10,x) = (46189 X^10-109395 X^8 + 90090 X^6 - 30030 X^4 + 3465 X^2-63)/256
194 
195  Recursion:
196 
197  P(0,x) = 1
198  P(1,x) = x
199  P(n,x) = ( (2*n-1)*x*P(n-1,x)-(n-1)*P(n-2,x) ) / n
200 
201  P'(0,x) = 0
202  P'(1,x) = 1
203  P'(N,x) = ( (2*N-1)*(P(N-1,x)+X*P'(N-1,x)-(N-1)*P'(N-2,x) ) / N
204 
205  Licensing:
206 
207  This code is distributed under the GNU LGPL license.
208 
209  Modified:
210 
211  08 August 2013
212 
213  Author:
214 
215  John Burkardt
216 
217  Reference:
218 
219  Milton Abramowitz, Irene Stegun,
220  Handbook of Mathematical Functions,
221  National Bureau of Standards, 1964,
222  ISBN: 0-486-61272-4,
223  LC: QA47.A34.
224 
225  Daniel Zwillinger, editor,
226  CRC Standard Mathematical Tables and Formulae,
227  30th Edition,
228  CRC Press, 1996.
229 
230  Parameters:
231 
232  Input, int M, the number of evaluation points.
233 
234  Input, int N, the highest order polynomial to evaluate.
235  Note that polynomials 0 through N will be evaluated.
236 
237  Input, double X[M], the evaluation points.
238 
239  Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre
240  polynomials of order 0 through N.
241  }}}*/
242 
243  int i,j;
244 
245  if(n<0) return NULL;
246 
247  IssmDouble* v = xNew<IssmDouble>(m*(n+1));
248 
249  for ( i = 0; i < m; i++ ) v[i*(n+1)+0] = 1.0;
250  if ( n < 1 ) return v;
251 
252  for ( i = 0; i < m; i++ ) v[i*(n+1)+1] = x[i];
253 
254  for ( i = 0; i < m; i++ ) {
255  for ( j = 2; j <= n; j++ ) {
256 
257  v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)]
258  - ( IssmDouble ) ( j - 1 ) * v[(j-2)+i*(n+1)] )
259  / ( IssmDouble ) ( j );
260 
261  }
262  }
263 
264  return v;
265 }
IssmDouble
double IssmDouble
Definition: types.h:37