Nektar++
Macros | Functions
Polylib_test.cpp File Reference
#include "Polylib.h"
#include <cmath>
#include <cstdio>

Go to the source code of this file.

Macros

#define NPLOWER   5
 
#define NPUPPER   15
 
#define EPS   1e-12
 
#define GAUSS_INT   1
 
#define GAUSS_RADAUM_INT   1
 
#define GAUSS_RADAUP_INT   1
 
#define GAUSS_LOBATTO_INT   1
 
#define GAUSS_DIFF   1
 
#define GAUSS_RADAUM_DIFF   1
 
#define GAUSS_RADAUP_DIFF   1
 
#define GAUSS_LOBATTO_DIFF   1
 
#define GAUSS_INTERP   1
 
#define GAUSS_RADAUM_INTERP   1
 
#define GAUSS_RADAUP_INTERP   1
 
#define GAUSS_LOBATTO_INTERP   1
 

Functions

double ddot (int, double *, int, double *, int)
 
double * dvector (int, int)
 
int test_gamma_fraction ()
 
int main ()
 

Macro Definition Documentation

◆ EPS

#define EPS   1e-12

Definition at line 130 of file Polylib_test.cpp.

◆ GAUSS_DIFF

#define GAUSS_DIFF   1

Definition at line 136 of file Polylib_test.cpp.

◆ GAUSS_INT

#define GAUSS_INT   1

Definition at line 132 of file Polylib_test.cpp.

◆ GAUSS_INTERP

#define GAUSS_INTERP   1

Definition at line 140 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_DIFF

#define GAUSS_LOBATTO_DIFF   1

Definition at line 139 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_INT

#define GAUSS_LOBATTO_INT   1

Definition at line 135 of file Polylib_test.cpp.

◆ GAUSS_LOBATTO_INTERP

#define GAUSS_LOBATTO_INTERP   1

Definition at line 143 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_DIFF

#define GAUSS_RADAUM_DIFF   1

Definition at line 137 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_INT

#define GAUSS_RADAUM_INT   1

Definition at line 133 of file Polylib_test.cpp.

◆ GAUSS_RADAUM_INTERP

#define GAUSS_RADAUM_INTERP   1

Definition at line 141 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_DIFF

#define GAUSS_RADAUP_DIFF   1

Definition at line 138 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_INT

#define GAUSS_RADAUP_INT   1

Definition at line 134 of file Polylib_test.cpp.

◆ GAUSS_RADAUP_INTERP

#define GAUSS_RADAUP_INTERP   1

Definition at line 142 of file Polylib_test.cpp.

◆ NPLOWER

#define NPLOWER   5

Definition at line 128 of file Polylib_test.cpp.

◆ NPUPPER

#define NPUPPER   15

Definition at line 129 of file Polylib_test.cpp.

Function Documentation

◆ ddot()

double ddot ( int  n,
double *  x,
int  incx,
double *  y,
int  incy 
)

Definition at line 738 of file Polylib_test.cpp.

739 {
740  register double sum = 0.;
741 
742  while (n--)
743  {
744  sum += (*x) * (*y);
745  x += incx;
746  y += incy;
747  }
748  return sum;
749 }

Referenced by main().

◆ dvector()

double * dvector ( int  nl,
int  nh 
)

Definition at line 751 of file Polylib_test.cpp.

752 {
753  double *v;
754 
755  v = (double *)malloc((unsigned)(nh - nl + 1) * sizeof(double));
756  return v - nl;
757 }

Referenced by main().

◆ main()

int main ( )

Definition at line 150 of file Polylib_test.cpp.

151 {
152  int np, n, i;
153  double *z, *w, *p, *q, sum = 0, alpha, beta, *d, *dt;
154 
155  z = dvector(0, NPUPPER - 1);
156  w = dvector(0, NPUPPER - 1);
157  p = dvector(0, NPUPPER - 1);
158 
159  q = dvector(0, NPUPPER * NPUPPER - 1);
160  d = dvector(0, NPUPPER * NPUPPER - 1);
161  dt = dvector(0, NPUPPER * NPUPPER - 1);
162 
163 #if GAUSS_INT
164  /* Gauss Integration */
165  printf("Begin checking Gauss integration\n");
166  alpha = -0.5;
167  while (alpha <= 5.0)
168  {
169  beta = -0.5;
170  while (beta <= 5.0)
171  {
172 
173  for (np = NPLOWER; np <= NPUPPER; ++np)
174  {
175  zwgj(z, w, np, alpha, beta);
176  for (n = 2; n < 2 * np - 1; ++n)
177  {
178  jacobfd(np, z, p, NULL, n, alpha, beta);
179  sum = ddot(np, w, 1, p, 1);
180  if (fabs(sum) > EPS)
181  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
182  "integal was %lg\n",
183  alpha, beta, np, n, sum);
184  }
185  }
186 
187  beta += 0.5;
188  }
189  printf("finished checking all beta values for alpha = %lf\n", alpha);
190  alpha += 0.5;
191  }
192  printf("Finished checking Gauss Integration\n");
193 #endif
194 
195 #if GAUSS_RADAUM_INT
196  /* Gauss Radau Integration */
197  printf("Begin checking Gauss Radau Integration\n");
198  alpha = -0.5;
199  while (alpha <= 5.0)
200  {
201  beta = -0.5;
202  while (beta <= 5.0)
203  {
204  for (np = NPLOWER; np <= NPUPPER; ++np)
205  {
206  zwgrjm(z, w, np, alpha, beta);
207  for (n = 2; n < 2 * np - 2; ++n)
208  {
209  jacobfd(np, z, p, NULL, n, alpha, beta);
210  sum = ddot(np, w, 1, p, 1);
211  if (fabs(sum) > EPS)
212  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
213  "integal was %lg\n",
214  alpha, beta, np, n, sum);
215  }
216  }
217 
218  beta += 0.5;
219  }
220  printf("finished checking all beta values for alpha = %lf\n", alpha);
221  alpha += 0.5;
222  }
223  printf("Finished checking Gauss Radau (z=-1) Integration\n");
224 #endif
225 
226 #if GAUSS_RADAUP_INT
227  /* Gauss Radau Integration */
228  printf("Begin checking Gauss Radau Integration\n");
229  alpha = -0.5;
230  while (alpha <= 5.0)
231  {
232  beta = -0.5;
233  while (beta <= 5.0)
234  {
235  for (np = NPLOWER; np <= NPUPPER; ++np)
236  {
237  zwgrjp(z, w, np, alpha, beta);
238  for (n = 2; n < 2 * np - 2; ++n)
239  {
240  jacobfd(np, z, p, NULL, n, alpha, beta);
241  sum = ddot(np, w, 1, p, 1);
242  if (fabs(sum) > EPS)
243  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
244  "integal was %lg\n",
245  alpha, beta, np, n, sum);
246  }
247  }
248 
249  beta += 0.5;
250  }
251  printf("finished checking all beta values for alpha = %lf\n", alpha);
252  alpha += 0.5;
253  }
254  printf("Finished checking Gauss Radau (z=1) Integration\n");
255 #endif
256 
257 #if GAUSS_LOBATTO_INT
258  /* Gauss Lobatto Integration */
259  printf("Begin checking Gauss Lobatto integration\n");
260  alpha = -0.5;
261  while (alpha <= 5.0)
262  {
263  beta = -0.5;
264  while (beta <= 5.0)
265  {
266 
267  for (np = NPLOWER; np <= NPUPPER; ++np)
268  {
269  zwglj(z, w, np, alpha, beta);
270  for (n = 2; n < 2 * np - 3; ++n)
271  {
272  jacobfd(np, z, p, NULL, n, alpha, beta);
273  sum = ddot(np, w, 1, p, 1);
274  if (fabs(sum) > EPS)
275  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
276  "integal was %lg\n",
277  alpha, beta, np, n, sum);
278  }
279  }
280 
281  beta += 0.5;
282  }
283  printf("finished checking all beta values for alpha = %lf\n", alpha);
284  alpha += 0.5;
285  }
286  printf("Finished checking Gauss Lobatto Integration\n");
287 #endif
288 
289 #if GAUSS_INT
290  printf("Begin checking integration through Gauss points\n");
291  alpha = -0.5;
292  while (alpha <= 5.0)
293  {
294  beta = -0.5;
295  while (beta <= 5.0)
296  {
297 
298  for (np = NPLOWER; np <= NPUPPER; ++np)
299  {
300  zwgj(z, w, np, alpha, beta);
301  Qg(q, z, np, 0);
302  for (n = 2; n < np - 1; ++n)
303  {
304  for (i = 0; i < np; ++i)
305  p[i] = pow(z[i], n);
306  sum = 0;
307  for (i = 0; i < np; ++i)
308  sum += fabs(ddot(np, q + i * np, 1, p, 1) -
309  pow(z[i], n + 1) / (n + 1));
310  sum /= np;
311  if (fabs(sum) > EPS)
312  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
313  "difference %lg\n",
314  alpha, beta, np, n, sum);
315  }
316  }
317  beta += 0.5;
318  }
319  printf("finished checking all beta values for alpha = %lf\n", alpha);
320  alpha += 0.5;
321  }
322  printf("Finished checking Gauss Jacobi integration\n");
323 #endif
324 
325 #if GAUSS_RADAUM_INT
326  printf("Begin checking integration through Gauss Radau points\n");
327  alpha = -0.5;
328  while (alpha <= 5.0)
329  {
330  beta = -0.5;
331  while (beta <= 5.0)
332  {
333 
334  for (np = NPLOWER; np <= NPUPPER; ++np)
335  {
336  zwgrjm(z, w, np, alpha, beta);
337  Qg(q, z, np, 0);
338  for (n = 2; n < np - 1; ++n)
339  {
340  for (i = 0; i < np; ++i)
341  p[i] = pow(z[i], n);
342  sum = 0;
343  for (i = 0; i < np; ++i)
344  sum += fabs(ddot(np, q + i * np, 1, p, 1) -
345  pow(z[i], n + 1) / (n + 1));
346  sum /= np;
347  if (fabs(sum) > EPS)
348  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
349  "difference %lg\n",
350  alpha, beta, np, n, sum);
351  }
352  }
353  beta += 0.5;
354  }
355  printf("finished checking all beta values for alpha = %lf\n", alpha);
356  alpha += 0.5;
357  }
358  printf("Finished checking Gauss Radau integration\n");
359 #endif
360 
361 #if GAUSS_RADAUP_INT
362  printf("Begin checking integration through Gauss Radau (z=1) points\n");
363  alpha = -0.5;
364  while (alpha <= 5.0)
365  {
366  beta = -0.5;
367  while (beta <= 5.0)
368  {
369 
370  for (np = NPLOWER; np <= NPUPPER; ++np)
371  {
372  zwgrjp(z, w, np, alpha, beta);
373  Qg(q, z, np, 0);
374  for (n = 2; n < np - 1; ++n)
375  {
376  for (i = 0; i < np; ++i)
377  p[i] = pow(z[i], n);
378  sum = 0;
379  for (i = 0; i < np; ++i)
380  sum += fabs(ddot(np, q + i * np, 1, p, 1) -
381  pow(z[i], n + 1) / (n + 1));
382  sum /= np;
383  if (fabs(sum) > EPS)
384  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
385  "difference %lg\n",
386  alpha, beta, np, n, sum);
387  }
388  }
389  beta += 0.5;
390  }
391  printf("finished checking all beta values for alpha = %lf\n", alpha);
392  alpha += 0.5;
393  }
394  printf("Finished checking Gauss Radau (z=1) integration\n");
395 #endif
396 
397 #if GAUSS_LOBATTO_INT
398  printf("Begin checking integration through Gauss Lobatto points\n");
399  alpha = -0.5;
400  while (alpha <= 5.0)
401  {
402  beta = -0.5;
403  while (beta <= 5.0)
404  {
405 
406  for (np = NPLOWER; np <= NPUPPER; ++np)
407  {
408  zwglj(z, w, np, alpha, beta);
409  Qg(q, z, np, 0);
410  for (n = 2; n < np - 1; ++n)
411  {
412  for (i = 0; i < np; ++i)
413  p[i] = pow(z[i], n);
414  sum = 0;
415  for (i = 0; i < np; ++i)
416  sum += fabs(ddot(np, q + i * np, 1, p, 1) -
417  pow(z[i], n + 1) / (n + 1));
418  sum /= np;
419  if (fabs(sum) > EPS)
420  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
421  "difference %lg\n",
422  alpha, beta, np, n, sum);
423  }
424  }
425  beta += 0.5;
426  }
427  printf("finished checking all beta values for alpha = %lf\n", alpha);
428  alpha += 0.5;
429  }
430  printf("Finished checking Gauss Lobatto integration\n");
431 #endif
432 
433 #if GAUSS_DIFF
434  printf("Begin checking differentiation through Gauss points\n");
435  alpha = -0.5;
436  while (alpha <= 5.0)
437  {
438  beta = -0.5;
439  while (beta <= 5.0)
440  {
441 
442  for (np = NPLOWER; np <= NPUPPER; ++np)
443  {
444  zwgj(z, w, np, alpha, beta);
445  Dgj(d, z, np, alpha, beta);
446  for (n = 2; n < np - 1; ++n)
447  {
448  for (i = 0; i < np; ++i)
449  p[i] = pow(z[i], n);
450  sum = 0;
451  for (i = 0; i < np; ++i)
452  sum += fabs(ddot(np, d + i, np, p, 1) -
453  n * pow(z[i], n - 1));
454  sum /= np;
455  if (fabs(sum) > EPS)
456  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
457  "difference %lg\n",
458  alpha, beta, np, n, sum);
459  }
460  }
461  beta += 0.5;
462  }
463  printf("finished checking all beta values for alpha = %lf\n", alpha);
464  alpha += 0.5;
465  }
466  printf("Finished checking Gauss Jacobi differentiation\n");
467 #endif
468 
469 #if GAUSS_RADAUM_DIFF
470  printf("Begin checking differentiation through Gauss Radau points\n");
471  alpha = -0.5;
472  while (alpha <= 5.0)
473  {
474  beta = -0.5;
475  while (beta <= 5.0)
476  {
477 
478  for (np = NPLOWER; np <= NPUPPER; ++np)
479  {
480  zwgrjm(z, w, np, alpha, beta);
481  Dgrjm(d, z, np, alpha, beta);
482  for (n = 2; n < np - 1; ++n)
483  {
484  for (i = 0; i < np; ++i)
485  p[i] = pow(z[i], n);
486  sum = 0;
487  for (i = 0; i < np; ++i)
488  sum += fabs(ddot(np, d + i, np, p, 1) -
489  n * pow(z[i], n - 1));
490  sum /= np;
491  if (fabs(sum) > EPS)
492  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
493  "difference %lg\n",
494  alpha, beta, np, n, sum);
495  }
496  }
497  beta += 0.5;
498  }
499  printf("finished checking all beta values for alpha = %lf\n", alpha);
500  alpha += 0.5;
501  }
502  printf("Finished checking Gauss Radau (z=-1) differentiation\n");
503 #endif
504 
505 #if GAUSS_RADAUP_DIFF
506  printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
507  alpha = -0.5;
508  while (alpha <= 5.0)
509  {
510  beta = -0.5;
511  while (beta <= 5.0)
512  {
513 
514  for (np = NPLOWER; np <= NPUPPER; ++np)
515  {
516  zwgrjp(z, w, np, alpha, beta);
517  Dgrjp(d, z, np, alpha, beta);
518  for (n = 2; n < np - 1; ++n)
519  {
520  for (i = 0; i < np; ++i)
521  p[i] = pow(z[i], n);
522  sum = 0;
523  for (i = 0; i < np; ++i)
524  sum += fabs(ddot(np, d + i, np, p, 1) -
525  n * pow(z[i], n - 1));
526  sum /= np;
527  if (fabs(sum) > EPS)
528  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
529  "difference %lg\n",
530  alpha, beta, np, n, sum);
531  }
532  }
533  beta += 0.5;
534  }
535  printf("finished checking all beta values for alpha = %lf\n", alpha);
536  alpha += 0.5;
537  }
538  printf("Finished checking Gauss Radau (z=1) differentiation\n");
539 #endif
540 
541 #if GAUSS_LOBATTO_DIFF
542  printf("Begin checking differentiation through Gauss Lobatto points\n");
543  alpha = -0.5;
544  while (alpha <= 5.0)
545  {
546  beta = -0.5;
547  while (beta <= 5.0)
548  {
549 
550  for (np = NPLOWER; np <= NPUPPER; ++np)
551  {
552  zwglj(z, w, np, alpha, beta);
553  Dglj(d, z, np, alpha, beta);
554  for (n = 2; n < np - 1; ++n)
555  {
556  for (i = 0; i < np; ++i)
557  p[i] = pow(z[i], n);
558  sum = 0;
559  for (i = 0; i < np; ++i)
560  sum += fabs(ddot(np, d + i, np, p, 1) -
561  n * pow(z[i], n - 1));
562  sum /= np;
563  if (fabs(sum) > EPS)
564  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
565  "difference %lg\n",
566  alpha, beta, np, n, sum);
567  }
568  }
569  beta += 0.5;
570  }
571  printf("finished checking all beta values for alpha = %lf\n", alpha);
572  alpha += 0.5;
573  }
574  printf("Finished checking Gauss Lobatto differentiation\n");
575 #endif
576 
577  /* check interpolation routines */
578 #if GAUSS_INTERP
579  printf("Begin checking interpolation through Gauss points\n");
580  alpha = -0.5;
581  while (alpha <= 5.0)
582  {
583  beta = -0.5;
584  while (beta <= 5.0)
585  {
586 
587  for (np = NPLOWER; np <= NPUPPER; ++np)
588  {
589  zwgj(z, w, np, alpha, beta);
590  for (n = 2; n < np - 1; ++n)
591  {
592  for (i = 0; i < np; ++i)
593  {
594  w[i] = 2.0 * i / (double)(np - 1) - 1.0;
595  p[i] = pow(z[i], n);
596  }
597  Imgj(d, z, w, np, np, alpha, beta);
598  sum = 0;
599  for (i = 0; i < np; ++i)
600  sum += fabs(ddot(np, d + i, np, p, 1) - pow(w[i], n));
601  sum /= np;
602  if (fabs(sum) > EPS)
603  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
604  "difference %lg\n",
605  alpha, beta, np, n, sum);
606  }
607  }
608  beta += 0.5;
609  }
610  printf("finished checking all beta values for alpha = %lf\n", alpha);
611  alpha += 0.5;
612  }
613  printf("Finished checking Gauss Jacobi interpolation\n");
614 #endif
615 
616 #if GAUSS_RADAUM_INTERP
617  printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
618  alpha = -0.5;
619  while (alpha <= 5.0)
620  {
621  beta = -0.5;
622  while (beta <= 5.0)
623  {
624 
625  for (np = NPLOWER; np <= NPUPPER; ++np)
626  {
627  zwgrjm(z, w, np, alpha, beta);
628  for (n = 2; n < np - 1; ++n)
629  {
630  for (i = 0; i < np; ++i)
631  {
632  w[i] = 2.0 * i / (double)(np - 1) - 1.0;
633  p[i] = pow(z[i], n);
634  }
635  Imgrjm(d, z, w, np, np, alpha, beta);
636  sum = 0;
637  for (i = 0; i < np; ++i)
638  sum += fabs(ddot(np, d + i, np, p, 1) - pow(w[i], n));
639  sum /= np;
640  if (fabs(sum) > EPS)
641  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
642  "difference %lg\n",
643  alpha, beta, np, n, sum);
644  }
645  }
646  beta += 0.5;
647  }
648  printf("finished checking all beta values for alpha = %lf\n", alpha);
649  alpha += 0.5;
650  }
651  printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
652 #endif
653 #if GAUSS_RADAUP_INTERP
654  printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
655  alpha = -0.5;
656  while (alpha <= 5.0)
657  {
658  beta = -0.5;
659  while (beta <= 5.0)
660  {
661 
662  for (np = NPLOWER; np <= NPUPPER; ++np)
663  {
664  zwgrjp(z, w, np, alpha, beta);
665  for (n = 2; n < np - 1; ++n)
666  {
667  for (i = 0; i < np; ++i)
668  {
669  w[i] = 2.0 * i / (double)(np - 1) - 1.0;
670  p[i] = pow(z[i], n);
671  }
672  Imgrjp(d, z, w, np, np, alpha, beta);
673  sum = 0;
674  for (i = 0; i < np; ++i)
675  sum += fabs(ddot(np, d + i, np, p, 1) - pow(w[i], n));
676  sum /= np;
677  if (fabs(sum) > EPS)
678  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
679  "difference %lg\n",
680  alpha, beta, np, n, sum);
681  }
682  }
683  beta += 0.5;
684  }
685  printf("finished checking all beta values for alpha = %lf\n", alpha);
686  alpha += 0.5;
687  }
688  printf("Finished checking Gauss Radau (z=1) interpolation\n");
689 #endif
690 
691 #if GAUSS_LOBATTO_INTERP
692  printf("Begin checking Interpolation through Gauss Lobatto points\n");
693  alpha = -0.5;
694  while (alpha <= 5.0)
695  {
696  beta = -0.5;
697  while (beta <= 5.0)
698  {
699 
700  for (np = NPLOWER; np <= NPUPPER; ++np)
701  {
702  zwglj(z, w, np, alpha, beta);
703  for (n = 2; n < np - 1; ++n)
704  {
705  for (i = 0; i < np; ++i)
706  {
707  w[i] = 2.0 * i / (double)(np - 1) - 1.0;
708  p[i] = pow(z[i], n);
709  }
710  Imglj(d, z, w, np, np, alpha, beta);
711  sum = 0;
712  for (i = 0; i < np; ++i)
713  sum += fabs(ddot(np, d + i, np, p, 1) - pow(w[i], n));
714  sum /= np;
715  if (fabs(sum) > EPS)
716  printf("alpha = %lf, beta = %lf, np = %d, n = %d "
717  "difference %lg\n",
718  alpha, beta, np, n, sum);
719  }
720  }
721  beta += 0.5;
722  }
723  printf("finished checking all beta values for alpha = %lf\n", alpha);
724  alpha += 0.5;
725  }
726  printf("Finished checking Gauss Lobatto interploation\n");
727 #endif
728  test_gamma_fraction();
729 
730  free(z);
731  free(w);
732  free(p);
733  free(d);
734  free(dt);
735  return 0;
736 }
dvector
double * dvector(int, int)
Definition: Polylib_test.cpp:751
NPUPPER
#define NPUPPER
Definition: Polylib_test.cpp:129
test_gamma_fraction
int test_gamma_fraction()
Definition: Polylib_test.cpp:759
EPS
#define EPS
Definition: Polylib_test.cpp:130
ddot
double ddot(int, double *, int, double *, int)
Definition: Polylib_test.cpp:738
NPLOWER
#define NPLOWER
Definition: Polylib_test.cpp:128
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Polylib::Dgj
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
Definition: Polylib.cpp:627
Polylib::zwgrjm
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
Definition: Polylib.cpp:161
Polylib::zwgrjp
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
Definition: Polylib.cpp:202
Polylib::Dglj
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
Definition: Polylib.cpp:787
Polylib::Imgrjp
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at ...
Definition: Polylib.cpp:1054
Polylib::Qg
void Qg(double *Q, const double *z, const int np, const int offset)
Compute the Integration Matrix.
Definition: Polylib.cpp:584
Polylib::zwglj
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
Definition: Polylib.cpp:241
Polylib::zwgj
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:133
Polylib::Imgj
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.
Definition: Polylib.cpp:992
Polylib::Imgrjm
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at...
Definition: Polylib.cpp:1023
Polylib::Dgrjm
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
Definition: Polylib.cpp:674
Polylib::Dgrjp
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.
Definition: Polylib.cpp:730
Polylib::Imglj
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.
Definition: Polylib.cpp:1085
Polylib::jacobfd
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1181

References Nektar::LibUtilities::beta, ddot(), Polylib::Dgj(), Polylib::Dglj(), Polylib::Dgrjm(), Polylib::Dgrjp(), dvector(), EPS, Polylib::Imgj(), Polylib::Imglj(), Polylib::Imgrjm(), Polylib::Imgrjp(), Polylib::jacobfd(), NPLOWER, NPUPPER, CellMLToNektar.cellml_metadata::p, Polylib::Qg(), test_gamma_fraction(), Polylib::zwgj(), Polylib::zwglj(), Polylib::zwgrjm(), and Polylib::zwgrjp().

◆ test_gamma_fraction()

int test_gamma_fraction ( )

Definition at line 759 of file Polylib_test.cpp.

760 {
761  double a = 362880.;
762  double c = Polylib::gammaFracGammaF(10, 0., 0, 1.);
763  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10., 1.,
764  c / a - 1.);
765  c = Polylib::gammaFracGammaF(1, 0., 12, -2.);
766  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 1., 10.,
767  c * a - 1.);
768 
769  a = 30.;
770  c = Polylib::gammaFracGammaF(4, 3., 5, 0.);
771  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5., 7.,
772  c / a - 1.);
773  c = Polylib::gammaFracGammaF(5, 0., 7, 0.);
774  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 7., 5.,
775  c * a - 1.);
776 
777  a = 21651.09375;
778  c = Polylib::gammaFracGammaF(7, 3.5, 5, 0.5);
779  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 5.5,
780  c / a - 1.);
781  c = Polylib::gammaFracGammaF(5, 0.5, 10, 0.5);
782  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5, 10.5,
783  c * a - 1.);
784 
785  a = 97429.921875;
786  c = Polylib::gammaFracGammaF(10, 0.5, 5, -0.5);
787  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 4.5,
788  c / a - 1.);
789  c = Polylib::gammaFracGammaF(5, -0.5, 11, -0.5);
790  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 4.5, 10.5,
791  c * a - 1.);
792 
793  a = 2279.0625;
794  c = Polylib::gammaFracGammaF(10, -0.5, 5, 0.5);
795  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 9.5, 5.5,
796  c / a - 1.);
797  c = Polylib::gammaFracGammaF(5, 0.5, 10, -0.5);
798  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5, 9.5,
799  c * a - 1.);
800 
801  a = 639383.8623046875;
802  c = Polylib::gammaFracGammaF(10, 0.5, 0, 0.5);
803  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 0.5,
804  c / a - 1.);
805  c = Polylib::gammaFracGammaF(0, 0.5, 10, 0.5);
806  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 0.5, 10.5,
807  c * a - 1.);
808 
809  a = 5967498288235982848.;
810  c = Polylib::gammaFracGammaF(100, 0., 90, 0.5);
811  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 100., 90.5,
812  c / a - 1.);
813  c = Polylib::gammaFracGammaF(90, 0.5, 100, 0.);
814  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 90.5, 100.,
815  c * a - 1.);
816 
817  a = 5618641603777298169856.;
818  c = Polylib::gammaFracGammaF(200, 0., 191, -0.5);
819  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200., 190.5,
820  c / a - 1.);
821  c = Polylib::gammaFracGammaF(190, 0.5, 200, 0.);
822  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.5, 200.,
823  c * a - 1.);
824 
825  a = 77396694214720029196288.;
826  c = Polylib::gammaFracGammaF(200, 0., 190, 0.);
827  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200., 190.,
828  c / a - 1.);
829  c = Polylib::gammaFracGammaF(190, 0., 200, 0.);
830  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190., 200.,
831  c * a - 1.);
832 
833  int Q = 2;
834  for (double alpha = -0.5; alpha <= 150.; alpha += 0.5)
835  {
836  for (double beta = -0.5; beta <= 150.; beta += 0.5)
837  {
838  a = Polylib::gammaF(Q + alpha) / Polylib::gammaF(Q + beta);
839  c = Polylib::gammaFracGammaF(Q, alpha, Q, beta) / a - 1.;
840  if (fabs(c) > 5.e-15)
841  {
842  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n",
843  Q + alpha, Q + beta, c);
844  }
845  }
846  }
847 
848  return 0;
849 }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1322
Polylib::gammaFracGammaF
double gammaFracGammaF(const int, const double, const int, const double)
Calculate fraction of two Gamma functions, , for integer values and halves.
Definition: Polylib.cpp:1371

References Nektar::LibUtilities::beta, Polylib::gammaF(), and Polylib::gammaFracGammaF().

Referenced by main().

Generated on Sat Mar 25 2023 12:35:54 for Nektar++ by doxygen 1.9.1