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