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