Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020 of Glen Hocky
3 :
4 : The FISST module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The FISST module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : //#ifndef __PLUMED_fisst_legendre_rule_fast_h
18 : //#define __PLUMED_fisst_legendre_rule_fast_h
19 : // adapted from:
20 : // https://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule_fast/legendre_rule_fast.html
21 : //
22 : #include "legendre_rule_fast.h"
23 :
24 : namespace PLMD {
25 : namespace fisst {
26 : //****************************************************************************80
27 :
28 2 : void legendre_compute_glr ( int n, double x[], double w[] )
29 :
30 : //****************************************************************************80
31 : //
32 : // Purpose:
33 : //
34 : // LEGENDRE_COMPUTE_GLR: Legendre quadrature by the Glaser-Liu-Rokhlin method.
35 : //
36 : // Licensing:
37 : //
38 : // This code is distributed under the GNU LGPL license.
39 : //
40 : // Modified:
41 : //
42 : // 20 October 2009
43 : //
44 : // Author:
45 : //
46 : // Original C++ version by Nick Hale.
47 : // This C++ version by John Burkardt.
48 : //
49 : // Reference:
50 : //
51 : // Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
52 : // A fast algorithm for the calculation of the roots of special functions,
53 : // SIAM Journal on Scientific Computing,
54 : // Volume 29, Number 4, pages 1420-1438, 2007.
55 : //
56 : // Parameters:
57 : //
58 : // Input, int N, the order.
59 : //
60 : // Output, double X[N], the abscissas.
61 : //
62 : // Output, double W[N], the weights.
63 : //
64 : {
65 : int i;
66 : double p;
67 : double pp;
68 : double w_sum;
69 : //
70 : // Get the value and derivative of the N-th Legendre polynomial at 0.
71 : //
72 2 : legendre_compute_glr0 ( n, &p, &pp );
73 : //
74 : // If N is odd, then zero is a root.
75 : //
76 2 : if ( n % 2 == 1 )
77 : {
78 2 : x[(n-1)/2] = p;
79 2 : w[(n-1)/2] = pp;
80 : }
81 : //
82 : // If N is even, we have to call a function to find the first root.
83 : //
84 : else
85 : {
86 0 : legendre_compute_glr2 ( p, n, &x[n/2], &w[n/2] );
87 : }
88 : //
89 : // Get the complete set of roots and derivatives.
90 : //
91 2 : legendre_compute_glr1 ( n, x, w );
92 : //
93 : // Compute the W.
94 : //
95 64 : for ( i = 0; i < n; i++ )
96 : {
97 62 : w[i] = 2.0 / ( 1.0 - x[i] ) / ( 1.0 + x[i] ) / w[i] / w[i];
98 : }
99 : w_sum = 0.0;
100 64 : for ( i = 0; i < n; i++ )
101 : {
102 62 : w_sum = w_sum + w[i];
103 : }
104 64 : for ( i = 0; i < n; i++ )
105 : {
106 62 : w[i] = 2.0 * w[i] / w_sum;
107 : }
108 2 : return;
109 : }
110 : //****************************************************************************80
111 :
112 2 : void legendre_compute_glr0 ( int n, double *p, double *pp )
113 :
114 : //****************************************************************************80
115 : //
116 : // Purpose:
117 : //
118 : // LEGENDRE_COMPUTE_GLR0 gets a starting value for the fast algorithm.
119 : //
120 : // Licensing:
121 : //
122 : // This code is distributed under the GNU LGPL license.
123 : //
124 : // Modified:
125 : //
126 : // 19 October 2009
127 : //
128 : // Author:
129 : //
130 : // Original C++ version by Nick Hale.
131 : // This C++ version by John Burkardt.
132 : //
133 : // Reference:
134 : //
135 : // Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
136 : // A fast algorithm for the calculation of the roots of special functions,
137 : // SIAM Journal on Scientific Computing,
138 : // Volume 29, Number 4, pages 1420-1438, 2007.
139 : //
140 : // Parameters:
141 : //
142 : // Input, int N, the order of the Legendre polynomial.
143 : //
144 : // Output, double *P, *PP, the value of the N-th Legendre polynomial
145 : // and its derivative at 0.
146 : //
147 : {
148 : double dk;
149 : int k;
150 : double pm1;
151 : double pm2;
152 : double ppm1;
153 : double ppm2;
154 :
155 : pm2 = 0.0;
156 : pm1 = 1.0;
157 : ppm2 = 0.0;
158 : ppm1 = 0.0;
159 :
160 64 : for ( k = 0; k < n; k++)
161 : {
162 62 : dk = ( double ) k;
163 62 : *p = - dk * pm2 / ( dk + 1.0 );
164 62 : *pp = ( ( 2.0 * dk + 1.0 ) * pm1 - dk * ppm2 ) / ( dk + 1.0 );
165 : pm2 = pm1;
166 62 : pm1 = *p;
167 : ppm2 = ppm1;
168 : ppm1 = *pp;
169 : }
170 2 : return;
171 : }
172 : //****************************************************************************80
173 :
174 2 : void legendre_compute_glr1 ( int n, double *x, double *w )
175 :
176 : //****************************************************************************80
177 : //
178 : // Purpose:
179 : //
180 : // LEGENDRE_COMPUTE_GLR1 gets the complete set of Legendre points and weights.
181 : //
182 : // Discussion:
183 : //
184 : // This routine requires that a starting estimate be provided for one
185 : // root and its derivative. This information will be stored in entry
186 : // (N+1)/2 if N is odd, or N/2 if N is even, of X and W.
187 : //
188 : // Licensing:
189 : //
190 : // This code is distributed under the GNU LGPL license.
191 : //
192 : // Modified:
193 : //
194 : // 19 October 2009
195 : //
196 : // Author:
197 : //
198 : // Original C++ version by Nick Hale.
199 : // This C++ version by John Burkardt.
200 : //
201 : // Reference:
202 : //
203 : // Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
204 : // A fast algorithm for the calculation of the roots of special functions,
205 : // SIAM Journal on Scientific Computing,
206 : // Volume 29, Number 4, pages 1420-1438, 2007.
207 : //
208 : // Parameters:
209 : //
210 : // Input, int N, the order of the Legendre polynomial.
211 : //
212 : // Input/output, double X[N]. On input, a starting value
213 : // has been set in one entry. On output, the roots of the Legendre
214 : // polynomial.
215 : //
216 : // Input/output, double W[N]. On input, a starting value
217 : // has been set in one entry. On output, the derivatives of the Legendre
218 : // polynomial at the zeros.
219 : //
220 : // Local Parameters:
221 : //
222 : // Local, int M, the number of terms in the Taylor expansion.
223 : //
224 : {
225 : double dk;
226 : double dn;
227 : double h;
228 : int j;
229 : int k;
230 : int l;
231 : int m = 30;
232 : int n2;
233 : static double pi = 3.141592653589793;
234 : int s;
235 : double *u;
236 : double *up;
237 : double xp;
238 :
239 2 : if ( n % 2 == 1 )
240 : {
241 2 : n2 = ( n - 1 ) / 2 - 1;
242 : s = 1;
243 : }
244 : else
245 : {
246 0 : n2 = n / 2 - 1;
247 : s = 0;
248 : }
249 :
250 2 : u = new double[m+2];
251 2 : up = new double[m+1];
252 :
253 2 : dn = ( double ) n;
254 :
255 32 : for ( j = n2 + 1; j < n - 1; j++ )
256 : {
257 30 : xp = x[j];
258 :
259 30 : h = rk2_leg ( pi/2.0, -pi/2.0, xp, n ) - xp;
260 :
261 30 : u[0] = 0.0;
262 30 : u[1] = 0.0;
263 30 : u[2] = w[j];
264 :
265 30 : up[0] = 0.0;
266 30 : up[1] = u[2];
267 :
268 900 : for ( k = 0; k <= m - 2; k++ )
269 : {
270 870 : dk = ( double ) k;
271 :
272 870 : u[k+3] =
273 : (
274 870 : 2.0 * xp * ( dk + 1.0 ) * u[k+2]
275 870 : + ( dk * ( dk + 1.0 ) - dn * ( dn + 1.0 ) ) * u[k+1] / ( dk + 1.0 )
276 870 : ) / ( 1.0 - xp ) / ( 1.0 + xp ) / ( dk + 2.0 );
277 :
278 870 : up[k+2] = ( dk + 2.0 ) * u[k+3];
279 : }
280 :
281 180 : for ( l = 0; l < 5; l++ )
282 : {
283 150 : h = h - ts_mult ( u, h, m ) / ts_mult ( up, h, m-1 );
284 : }
285 :
286 30 : x[j+1] = xp + h;
287 30 : w[j+1] = ts_mult ( up, h, m - 1 );
288 : }
289 :
290 34 : for ( k = 0; k <= n2 + s; k++ )
291 : {
292 32 : x[k] = - x[n-1-k];
293 32 : w[k] = w[n-1-k];
294 : }
295 2 : return;
296 : }
297 : //****************************************************************************80
298 :
299 0 : void legendre_compute_glr2 ( double pn0, int n, double *x1, double *d1 )
300 :
301 : //****************************************************************************80
302 : //
303 : // Purpose:
304 : //
305 : // LEGENDRE_COMPUTE_GLR2 finds the first real root.
306 : //
307 : // Discussion:
308 : //
309 : // This function is only called if N is even.
310 : //
311 : // Licensing:
312 : //
313 : // This code is distributed under the GNU LGPL license.
314 : //
315 : // Modified:
316 : //
317 : // 19 October 2009
318 : //
319 : // Author:
320 : //
321 : // Original C++ version by Nick Hale.
322 : // This C++ version by John Burkardt.
323 : //
324 : // Reference:
325 : //
326 : // Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
327 : // A fast algorithm for the calculation of the roots of special functions,
328 : // SIAM Journal on Scientific Computing,
329 : // Volume 29, Number 4, pages 1420-1438, 2007.
330 : //
331 : // Parameters:
332 : //
333 : // Input, double PN0, the value of the N-th Legendre polynomial
334 : // at 0.
335 : //
336 : // Input, int N, the order of the Legendre polynomial.
337 : //
338 : // Output, double *X1, the first real root.
339 : //
340 : // Output, double *D1, the derivative at X1.
341 : //
342 : // Local Parameters:
343 : //
344 : // Local, int M, the number of terms in the Taylor expansion.
345 : //
346 : {
347 : double dk;
348 : double dn;
349 : int k;
350 : int l;
351 : int m = 30;
352 : static double pi = 3.141592653589793;
353 : double t;
354 : double *u;
355 : double *up;
356 :
357 : t = 0.0;
358 0 : *x1 = rk2_leg ( t, -pi/2.0, 0.0, n );
359 :
360 0 : u = new double[m+2];
361 0 : up = new double[m+1];
362 :
363 0 : dn = ( double ) n;
364 : //
365 : // U[0] and UP[0] are never used.
366 : // U[M+1] is set, but not used, and UP[M] is set and not used.
367 : // What gives?
368 : //
369 0 : u[0] = 0.0;
370 0 : u[1] = pn0;
371 :
372 0 : up[0] = 0.0;
373 :
374 0 : for ( k = 0; k <= m - 2; k = k + 2 )
375 : {
376 0 : dk = ( double ) k;
377 :
378 0 : u[k+2] = 0.0;
379 0 : u[k+3] = ( dk * ( dk + 1.0 ) - dn * ( dn + 1.0 ) ) * u[k+1]
380 0 : / (dk + 1.0) / (dk + 2.0 );
381 :
382 0 : up[k+1] = 0.0;
383 0 : up[k+2] = ( dk + 2.0 ) * u[k+3];
384 : }
385 :
386 0 : for ( l = 0; l < 5; l++ )
387 : {
388 0 : *x1 = *x1 - ts_mult ( u, *x1, m ) / ts_mult ( up, *x1, m-1 );
389 : }
390 0 : *d1 = ts_mult ( up, *x1, m-1 );
391 :
392 0 : return;
393 : }
394 :
395 : //****************************************************************************80
396 :
397 2 : void rescale ( double a, double b, int n, double x[], double w[] )
398 :
399 : //****************************************************************************80
400 : //
401 : // Purpose:
402 : //
403 : // RESCALE rescales a Legendre quadrature rule from [-1,+1] to [A,B].
404 : //
405 : // Licensing:
406 : //
407 : // This code is distributed under the GNU LGPL license.
408 : //
409 : // Modified:
410 : //
411 : // 18 October 2009
412 : //
413 : // Author:
414 : //
415 : // Original MATLAB version by Nick Hale.
416 : // C++ version by John Burkardt.
417 : //
418 : // Reference:
419 : //
420 : // Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
421 : // A fast algorithm for the calculation of the roots of special functions,
422 : // SIAM Journal on Scientific Computing,
423 : // Volume 29, Number 4, pages 1420-1438, 2007.
424 : //
425 : // Parameters:
426 : //
427 : // Input, double A, B, the endpoints of the new interval.
428 : //
429 : // Input, int N, the order.
430 : //
431 : // Input/output, double X[N], on input, the abscissas for [-1,+1].
432 : // On output, the abscissas for [A,B].
433 : //
434 : // Input/output, double W[N], on input, the weights for [-1,+1].
435 : // On output, the weights for [A,B].
436 : //
437 : {
438 : int i;
439 :
440 64 : for ( i = 0; i < n; i++ )
441 : {
442 62 : x[i] = ( ( a + b ) + ( b - a ) * x[i] ) / 2.0;
443 : }
444 64 : for ( i = 0; i < n; i++ )
445 : {
446 62 : w[i] = ( b - a ) * w[i] / 2.0;
447 : }
448 2 : return;
449 : }
450 : //****************************************************************************80
451 :
452 30 : double rk2_leg ( double t1, double t2, double x, int n )
453 :
454 : //****************************************************************************80
455 : //
456 : // Purpose:
457 : //
458 : // RK2_LEG advances the value of X(T) using a Runge-Kutta method.
459 : //
460 : // Licensing:
461 : //
462 : // This code is distributed under the GNU LGPL license.
463 : //
464 : // Modified:
465 : //
466 : // 22 October 2009
467 : //
468 : // Author:
469 : //
470 : // Original C++ version by Nick Hale.
471 : // This C++ version by John Burkardt.
472 : //
473 : // Parameters:
474 : //
475 : // Input, double T1, T2, the range of the integration interval.
476 : //
477 : // Input, double X, the value of X at T1.
478 : //
479 : // Input, int N, the number of steps to take.
480 : //
481 : // Output, double RK2_LEG, the value of X at T2.
482 : //
483 : {
484 : double f;
485 : double h;
486 : int j;
487 : double k1;
488 : double k2;
489 : int m = 10;
490 : double snn1;
491 : double t;
492 :
493 30 : h = ( t2 - t1 ) / ( double ) m;
494 30 : snn1 = sqrt ( ( double ) ( n * ( n + 1 ) ) );
495 : t = t1;
496 :
497 330 : for ( j = 0; j < m; j++ )
498 : {
499 300 : f = ( 1.0 - x ) * ( 1.0 + x );
500 300 : k1 = - h * f / ( snn1 * sqrt ( f ) - 0.5 * x * sin ( 2.0 * t ) );
501 300 : x = x + k1;
502 :
503 300 : t = t + h;
504 :
505 300 : f = ( 1.0 - x ) * ( 1.0 + x );
506 300 : k2 = - h * f / ( snn1 * sqrt ( f ) - 0.5 * x * sin ( 2.0 * t ) );
507 300 : x = x + 0.5 * ( k2 - k1 );
508 : }
509 30 : return x;
510 : }
511 : //****************************************************************************80
512 :
513 330 : double ts_mult ( double *u, double h, int n )
514 :
515 : //****************************************************************************80
516 : //
517 : // Purpose:
518 : //
519 : // TS_MULT evaluates a polynomial.
520 : //
521 : // Licensing:
522 : //
523 : // This code is distributed under the GNU LGPL license.
524 : //
525 : // Modified:
526 : //
527 : // 17 May 2013
528 : //
529 : // Author:
530 : //
531 : // Original C++ version by Nick Hale.
532 : // This C++ version by John Burkardt.
533 : //
534 : // Parameters:
535 : //
536 : // Input, double U[N+1], the polynomial coefficients.
537 : // U[0] is ignored.
538 : //
539 : // Input, double H, the polynomial argument.
540 : //
541 : // Input, int N, the number of terms to compute.
542 : //
543 : // Output, double TS_MULT, the value of the polynomial.
544 : //
545 : {
546 : double hk;
547 : int k;
548 : double ts;
549 :
550 : ts = 0.0;
551 : hk = 1.0;
552 10050 : for ( k = 1; k<= n; k++ )
553 : {
554 9720 : ts = ts + u[k] * hk;
555 9720 : hk = hk * h;
556 : }
557 330 : return ts;
558 : }
559 :
560 : //close the namespaces
561 : }
562 : }
563 :
564 : //#endif
565 :
|
