Nektar++
Polylib_test.cpp
Go to the documentation of this file.
1 #include "Polylib.h"
2 #include <cstdio>
3 #include <cmath>
4 
5 using namespace std;
6 using namespace Polylib;
7 
8 /* --------------------------------------------------------------------
9  To compile:
10  g++ -g -c Polylib.cpp -I ./ -I ../../
11  g++ -g -c Polylib_test.cpp -I ./ -I ../../
12  g++ -g -o polytest Polylib_test.o Polylib.o -lm
13  * --------------------------------------------------------------------*/
14 
15 /* -------------------------------------------------------------------
16  This is a routine to test the integration, differentiation and
17  interpolation routines in the polylib.c.
18 
19  First, it performs the integral
20 
21  /1 alpha beta alpha,beta
22  | (1-x) (1+x) P (x) dx = 0
23  /-1 n
24 
25  for all -0.5 <= alpha <= 5 (increments of 0.5)
26  -0.5 <= beta <= 5 (increments of 0.5)
27 
28  using np points where
29  NPLOWER <= np <= NPUPPER
30  2 <= n <= 2*np - delta
31 
32  delta = 1 (gauss), 2(radau), 3(lobatto).
33  The integral is evaluated and if it is larger that EPS then the
34  value of alpha,beta,np,n and the integral is printed to the screen.
35 
36  After every alpha value the statement
37  "finished checking all beta values for alpha = #"
38  is printed
39 
40  The routine then evaluates the derivate of
41 
42  d n n-1
43  -- x = n x
44  dx
45 
46  for all -0.5 <= alpha <= 5 (increments of 0.5)
47  -0.5 <= beta <= 5 (increments of 0.5)
48 
49  using np points where
50  NPLOWER <= np <= NPUPPER
51  2 <= n <= np - 1
52 
53  The error is check in a pointwise sense and if it is larger than
54  EPS then the value of alpha,beta,np,n and the error is printed to
55  the screen. After every alpha value the statement
56  "finished checking all beta values for alpha = #"
57  is printed
58 
59  Finally the routine evaluates the interpolation of
60 
61  n n
62  z to x
63 
64  where z are the quadrature zeros and x are the equispaced points
65 
66  2*i
67  x = ----- - 1.0 (0 <= i <= np-1)
68  i (np-1)
69 
70 
71  for all -0.5 <= alpha <= 5 (increments of 0.5)
72  -0.5 <= beta <= 5 (increments of 0.5)
73 
74  using np points where
75  NPLOWER <= np <= NPUPPER
76  2 <= n <= np - 1
77 
78  The error is check in a pointwise sense and if it is larger than
79  EPS then the value of alpha,beta,np,n and the error is printed to
80  the screen. After every alpha value the statement
81  "finished checking all beta values for alpha = #"
82  is printed
83 
84  The above checks are performed for all the Gauss, Gauss-Radau and
85  Gauss-Lobatto points. If you want to disable any routine then set
86  GAUSS_INT, GAUSS_RADAU_INT, GAUSS_LOBATTO_INT = 0
87  for the integration rouintes
88  GAUSS_DIFF,GAUSS_RADAU_DIFF, GAUSS_LOBATTO_DIFF = 0
89  for the differentiation routines
90  GAUSS_INTERP,GAUSS_RADAU_INTERP, GAUSS_LOBATTO_INTERP = 0
91  for the interpolation routines.
92  ------------------------------------------------------------------*/
93 
94 #define NPLOWER 5
95 #define NPUPPER 15
96 #define EPS 1e-12
97 
98 #define GAUSS_INT 1
99 #define GAUSS_RADAUM_INT 1
100 #define GAUSS_RADAUP_INT 1
101 #define GAUSS_LOBATTO_INT 1
102 #define GAUSS_DIFF 1
103 #define GAUSS_RADAUM_DIFF 1
104 #define GAUSS_RADAUP_DIFF 1
105 #define GAUSS_LOBATTO_DIFF 1
106 #define GAUSS_INTERP 1
107 #define GAUSS_RADAUM_INTERP 1
108 #define GAUSS_RADAUP_INTERP 1
109 #define GAUSS_LOBATTO_INTERP 1
110 
111 /* local routines */
112 double ddot (int, double *, int, double *, int);
113 double *dvector (int, int);
114 int test_gamma_fraction();
115 
116 int main(){
117  int np,n,i;
118  double *z,*w,*p,sum=0,alpha,beta,*d,*dt;
119 
120  z = dvector(0,NPUPPER-1);
121  w = dvector(0,NPUPPER-1);
122  p = dvector(0,NPUPPER-1);
123 
124  d = dvector(0,NPUPPER*NPUPPER-1);
125  dt = dvector(0,NPUPPER*NPUPPER-1);
126 
127 #if GAUSS_INT
128  /* Gauss Integration */
129  printf("Begin checking Gauss integration\n");
130  alpha = -0.5;
131  while(alpha <= 5.0){
132  beta = -0.5;
133  while(beta <= 5.0){
134 
135  for(np = NPLOWER; np <= NPUPPER; ++np){
136  zwgj(z,w,np,alpha,beta);
137  for(n = 2; n < 2*np-1; ++n){
138  jacobfd(np,z,p,NULL,n,alpha,beta);
139  sum = ddot(np,w,1,p,1);
140  if(fabs(sum)>EPS)
141  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
142  ,alpha,beta,np,n,sum);
143  }
144  }
145 
146  beta += 0.5;
147  }
148  printf("finished checking all beta values for alpha = %lf\n",alpha);
149  alpha += 0.5;
150  }
151  printf("Finished checking Gauss Integration\n");
152 #endif
153 
154 #if GAUSS_RADAUM_INT
155  /* Gauss Radau Integration */
156  printf("Begin checking Gauss Radau Integration\n");
157  alpha = -0.5;
158  while(alpha <= 5.0){
159  beta = -0.5;
160  while(beta <= 5.0){
161  for(np = NPLOWER; np <= NPUPPER; ++np){
162  zwgrjm(z,w,np,alpha,beta);
163  for(n = 2; n < 2*np-2; ++n){
164  jacobfd(np,z,p,NULL,n,alpha,beta);
165  sum = ddot(np,w,1,p,1);
166  if(fabs(sum)>EPS)
167  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
168  ,alpha,beta,np,n,sum);
169  }
170  }
171 
172  beta += 0.5;
173  }
174  printf("finished checking all beta values for alpha = %lf\n",alpha);
175  alpha += 0.5;
176  }
177  printf("Finished checking Gauss Radau (z=-1) Integration\n");
178 #endif
179 
180 
181 #if GAUSS_RADAUP_INT
182  /* Gauss Radau Integration */
183  printf("Begin checking Gauss Radau Integration\n");
184  alpha = -0.5;
185  while(alpha <= 5.0){
186  beta = -0.5;
187  while(beta <= 5.0){
188  for(np = NPLOWER; np <= NPUPPER; ++np){
189  zwgrjp(z,w,np,alpha,beta);
190  for(n = 2; n < 2*np-2; ++n){
191  jacobfd(np,z,p,NULL,n,alpha,beta);
192  sum = ddot(np,w,1,p,1);
193  if(fabs(sum)>EPS)
194  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
195  ,alpha,beta,np,n,sum);
196  }
197  }
198 
199  beta += 0.5;
200  }
201  printf("finished checking all beta values for alpha = %lf\n",alpha);
202  alpha += 0.5;
203  }
204  printf("Finished checking Gauss Radau (z=1) Integration\n");
205 #endif
206 
207 #if GAUSS_LOBATTO_INT
208  /* Gauss Lobatto Integration */
209  printf("Begin checking Gauss Lobatto integration\n");
210  alpha = -0.5;
211  while(alpha <= 5.0){
212  beta = -0.5;
213  while(beta <= 5.0){
214 
215  for(np = NPLOWER; np <= NPUPPER; ++np){
216  zwglj(z,w,np,alpha,beta);
217  for(n = 2; n < 2*np-3; ++n){
218  jacobfd(np,z,p,NULL,n,alpha,beta);
219  sum = ddot(np,w,1,p,1);
220  if(fabs(sum)>EPS)
221  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
222  ,alpha,beta,np,n,sum);
223  }
224  }
225 
226  beta += 0.5;
227  }
228  printf("finished checking all beta values for alpha = %lf\n",alpha);
229  alpha += 0.5;
230  }
231  printf("Finished checking Gauss Lobatto Integration\n");
232 #endif
233 
234 #if GAUSS_DIFF
235  printf("Begin checking differentiation through Gauss points\n");
236  alpha = -0.5;
237  while(alpha <= 5.0){
238  beta = -0.5;
239  while(beta <= 5.0){
240 
241  for(np = NPLOWER; np <= NPUPPER; ++np){
242  zwgj(z,w,np,alpha,beta);
243  for(n = 2; n < np-1; ++n){
244  Dgj(d,z,np,alpha,beta);
245 
246  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
247  sum = 0;
248  for(i = 0; i < np; ++i)
249  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
250  sum /= np;
251  if(fabs(sum)>EPS)
252  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
253  ,alpha,beta,np,n,sum);
254  }
255  }
256  beta += 0.5;
257  }
258  printf("finished checking all beta values for alpha = %lf\n",alpha);
259  alpha += 0.5;
260  }
261  printf("Finished checking Gauss Jacobi differentiation\n");
262 #endif
263 
264 #if GAUSS_RADAUM_DIFF
265  printf("Begin checking differentiation through Gauss Radau points\n");
266  alpha = -0.5;
267  while(alpha <= 5.0){
268  beta = -0.5;
269  while(beta <= 5.0){
270 
271  for(np = NPLOWER; np <= NPUPPER; ++np){
272  zwgrjm(z,w,np,alpha,beta);
273  for(n = 2; n < np-1; ++n){
274  Dgrjm(d,z,np,alpha,beta);
275  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
276  sum = 0;
277  for(i = 0; i < np; ++i)
278  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
279  sum /= np;
280  if(fabs(sum)>EPS)
281  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
282  ,alpha,beta,np,n,sum);
283  }
284  }
285  beta += 0.5;
286  }
287  printf("finished checking all beta values for alpha = %lf\n",alpha);
288  alpha += 0.5;
289  }
290  printf("Finished checking Gauss Radau (z=-1) differentiation\n");
291 #endif
292 
293 #if GAUSS_RADAUP_DIFF
294  printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
295  alpha = -0.5;
296  while(alpha <= 5.0){
297  beta = -0.5;
298  while(beta <= 5.0){
299 
300  for(np = NPLOWER; np <= NPUPPER; ++np){
301  zwgrjp(z,w,np,alpha,beta);
302  for(n = 2; n < np-1; ++n){
303  Dgrjp(d,z,np,alpha,beta);
304  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
305  sum = 0;
306  for(i = 0; i < np; ++i)
307  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
308  sum /= np;
309  if(fabs(sum)>EPS)
310  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
311  ,alpha,beta,np,n,sum);
312  }
313  }
314  beta += 0.5;
315  }
316  printf("finished checking all beta values for alpha = %lf\n",alpha);
317  alpha += 0.5;
318  }
319  printf("Finished checking Gauss Radau (z=1) differentiation\n");
320 #endif
321 
322 #if GAUSS_LOBATTO_DIFF
323  printf("Begin checking differentiation through Gauss Lobatto points\n");
324  alpha = -0.5;
325  while(alpha <= 5.0){
326  beta = -0.5;
327  while(beta <= 5.0){
328 
329  for(np = NPLOWER; np <= NPUPPER; ++np){
330  zwglj(z,w,np,alpha,beta);
331  for(n = 2; n < np-1; ++n){
332  Dglj(d,z,np,alpha,beta);
333  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
334  sum = 0;
335  for(i = 0; i < np; ++i)
336  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
337  sum /= np;
338  if(fabs(sum)>EPS)
339  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
340  ,alpha,beta,np,n,sum);
341  }
342  }
343  beta += 0.5;
344  }
345  printf("finished checking all beta values for alpha = %lf\n",alpha);
346  alpha += 0.5;
347  }
348  printf("Finished checking Gauss Lobatto differentiation\n");
349 #endif
350 
351  /* check interpolation routines */
352 #if GAUSS_INTERP
353  printf("Begin checking interpolation through Gauss points\n");
354  alpha = -0.5;
355  while(alpha <= 5.0){
356  beta = -0.5;
357  while(beta <= 5.0){
358 
359  for(np = NPLOWER; np <= NPUPPER; ++np){
360  zwgj(z,w,np,alpha,beta);
361  for(n = 2; n < np-1; ++n){
362  for(i = 0; i < np; ++i) {
363  w[i] = 2.0*i/(double)(np-1)-1.0;
364  p[i] = pow(z[i],n);
365  }
366  Imgj(d,z,w,np,np,alpha,beta);
367  sum = 0;
368  for(i = 0; i < np; ++i)
369  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
370  sum /= np;
371  if(fabs(sum)>EPS)
372  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
373  ,alpha,beta,np,n,sum);
374  }
375  }
376  beta += 0.5;
377  }
378  printf("finished checking all beta values for alpha = %lf\n",alpha);
379  alpha += 0.5;
380  }
381  printf("Finished checking Gauss Jacobi interpolation\n");
382 #endif
383 
384 #if GAUSS_RADAUM_INTERP
385  printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
386  alpha = -0.5;
387  while(alpha <= 5.0){
388  beta = -0.5;
389  while(beta <= 5.0){
390 
391  for(np = NPLOWER; np <= NPUPPER; ++np){
392  zwgrjm(z,w,np,alpha,beta);
393  for(n = 2; n < np-1; ++n){
394  for(i = 0; i < np; ++i) {
395  w[i] = 2.0*i/(double)(np-1)-1.0;
396  p[i] = pow(z[i],n);
397  }
398  Imgrjm(d,z,w,np,np,alpha,beta);
399  sum = 0;
400  for(i = 0; i < np; ++i)
401  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
402  sum /= np;
403  if(fabs(sum)>EPS)
404  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
405  ,alpha,beta,np,n,sum);
406  }
407  }
408  beta += 0.5;
409  }
410  printf("finished checking all beta values for alpha = %lf\n",alpha);
411  alpha += 0.5;
412  }
413  printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
414 #endif
415 #if GAUSS_RADAUP_INTERP
416  printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
417  alpha = -0.5;
418  while(alpha <= 5.0){
419  beta = -0.5;
420  while(beta <= 5.0){
421 
422  for(np = NPLOWER; np <= NPUPPER; ++np){
423  zwgrjp(z,w,np,alpha,beta);
424  for(n = 2; n < np-1; ++n){
425  for(i = 0; i < np; ++i) {
426  w[i] = 2.0*i/(double)(np-1)-1.0;
427  p[i] = pow(z[i],n);
428  }
429  Imgrjp(d,z,w,np,np,alpha,beta);
430  sum = 0;
431  for(i = 0; i < np; ++i)
432  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
433  sum /= np;
434  if(fabs(sum)>EPS)
435  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
436  ,alpha,beta,np,n,sum);
437  }
438  }
439  beta += 0.5;
440  }
441  printf("finished checking all beta values for alpha = %lf\n",alpha);
442  alpha += 0.5;
443  }
444  printf("Finished checking Gauss Radau (z=1) interpolation\n");
445 #endif
446 
447 #if GAUSS_LOBATTO_INTERP
448  printf("Begin checking Interpolation through Gauss Lobatto points\n");
449  alpha = -0.5;
450  while(alpha <= 5.0){
451  beta = -0.5;
452  while(beta <= 5.0){
453 
454  for(np = NPLOWER; np <= NPUPPER; ++np){
455  zwglj(z,w,np,alpha,beta);
456  for(n = 2; n < np-1; ++n){
457  for(i = 0; i < np; ++i) {
458  w[i] = 2.0*i/(double)(np-1)-1.0;
459  p[i] = pow(z[i],n);
460  }
461  Imglj(d,z,w,np,np,alpha,beta);
462  sum = 0;
463  for(i = 0; i < np; ++i)
464  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
465  sum /= np;
466  if(fabs(sum)>EPS)
467  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
468  ,alpha,beta,np,n,sum);
469  }
470  }
471  beta += 0.5;
472  }
473  printf("finished checking all beta values for alpha = %lf\n",alpha);
474  alpha += 0.5;
475  }
476  printf("Finished checking Gauss Lobatto interploation\n");
477 #endif
478  test_gamma_fraction();
479 
480 
481  free(z); free(w); free(p); free(d); free(dt);
482  return 0;
483 }
484 
485 double ddot (int n, double *x, int incx, double *y, int incy)
486 {
487  register double sum = 0.;
488 
489  while (n--) {
490  sum += (*x) * (*y);
491  x += incx;
492  y += incy;
493  }
494  return sum;
495 }
496 
497 
498 double *dvector(int nl,int nh)
499 {
500  double *v;
501 
502  v = (double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
503  return v-nl;
504 }
505 
506 int test_gamma_fraction(){
507  double a = 362880.;
508  double c = Polylib::gammaFracGammaF(10, 0., 0, 1.);
509  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.,1., c/a-1.);
510  c = Polylib::gammaFracGammaF(1, 0., 12, -2.);
511  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 1.,10., c*a-1.);
512 
513  a = 30.;
514  c = Polylib::gammaFracGammaF(4, 3., 5, 0.);
515  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.,7., c/a-1.);
516  c = Polylib::gammaFracGammaF(5, 0., 7, 0.);
517  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 7.,5., c*a-1.);
518 
519  a = 21651.09375;
520  c = Polylib::gammaFracGammaF(7, 3.5, 5, 0.5);
521  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,5.5, c/a-1.);
522  c = Polylib::gammaFracGammaF(5, 0.5, 10, 0.5);
523  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5,10.5, c*a-1.);
524 
525  a = 97429.921875;
526  c = Polylib::gammaFracGammaF(10, 0.5, 5, -0.5);
527  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,4.5, c/a-1.);
528  c = Polylib::gammaFracGammaF(5, -0.5, 11, -0.5);
529  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 4.5,10.5, c*a-1.);
530 
531  a = 2279.0625;
532  c = Polylib::gammaFracGammaF(10, -0.5, 5, 0.5);
533  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 9.5,5.5, c/a-1.);
534  c = Polylib::gammaFracGammaF(5, 0.5, 10, -0.5);
535  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5,9.5, c*a-1.);
536 
537  a = 639383.8623046875;
538  c = Polylib::gammaFracGammaF(10, 0.5, 0, 0.5);
539  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5,0.5, c/a-1.);
540  c = Polylib::gammaFracGammaF(0, 0.5, 10, 0.5);
541  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 0.5,10.5, c*a-1.);
542 
543  a = 5967498288235982848.;
544  c = Polylib::gammaFracGammaF(100, 0., 90, 0.5);
545  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 100.,90.5, c/a-1.);
546  c = Polylib::gammaFracGammaF(90, 0.5, 100, 0.);
547  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 90.5,100., c*a-1.);
548 
549  a = 5618641603777298169856.;
550  c = Polylib::gammaFracGammaF(200, 0., 191, -0.5);
551  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200.,190.5, c/a-1.);
552  c = Polylib::gammaFracGammaF(190, 0.5, 200, 0.);
553  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.5,200., c*a-1.);
554 
555  a = 77396694214720029196288.;
556  c = Polylib::gammaFracGammaF(200, 0., 190, 0.);
557  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200.,190., c/a-1.);
558  c = Polylib::gammaFracGammaF(190, 0., 200, 0.);
559  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.,200., c*a-1.);
560 
561  int Q=2;
562  for(double alpha=-0.5; alpha<=150.; alpha+=0.5)
563  {
564  for(double beta=-0.5; beta<=150.; beta+=0.5)
565  {
566  a = Polylib::gammaF(Q+alpha)/Polylib::gammaF(Q+beta);
567  c = Polylib::gammaFracGammaF(Q, alpha, Q, beta) / a - 1.;
568  if(fabs(c)>5.e-15)
569  {
570  printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", Q+alpha,Q+beta, c);
571  }
572  }
573  }
574 
575  return 0;
576 }
dvector
double * dvector(int, int)
Definition: Polylib_test.cpp:498
NPUPPER
#define NPUPPER
Definition: Polylib_test.cpp:95
test_gamma_fraction
int test_gamma_fraction()
Definition: Polylib_test.cpp:506
EPS
#define EPS
Definition: Polylib_test.cpp:96
ddot
double ddot(int, double *, int, double *, int)
Definition: Polylib_test.cpp:485
NPLOWER
#define NPLOWER
Definition: Polylib_test.cpp:94
main
int main()
Definition: Polylib_test.cpp:116
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Polylib
The namespace associated with the the Polylib library (Polylib introduction)
Definition: Polylib.cpp:17
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1161
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:546
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:119
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:158
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:690
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:946
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:195
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:91
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:887
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:1203
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:916
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:589
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:639
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:977
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:1034
Generated on Wed Nov 24 2021 10:33:28 for Nektar++ by doxygen 1.9.1