121 double *z,*w,*
p,sum=0,alpha,beta,*d,*dt;
132 printf(
"Begin checking Gauss integration\n");
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);
144 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
145 ,alpha,beta,np,n,sum);
151 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
154 printf(
"Finished checking Gauss Integration\n");
159 printf(
"Begin checking Gauss Radau Integration\n");
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);
170 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
171 ,alpha,beta,np,n,sum);
177 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
180 printf(
"Finished checking Gauss Radau (z=-1) Integration\n");
186 printf(
"Begin checking Gauss Radau Integration\n");
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);
197 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
198 ,alpha,beta,np,n,sum);
204 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
207 printf(
"Finished checking Gauss Radau (z=1) Integration\n");
210 #if GAUSS_LOBATTO_INT
212 printf(
"Begin checking Gauss Lobatto integration\n");
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);
224 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
225 ,alpha,beta,np,n,sum);
231 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
234 printf(
"Finished checking Gauss Lobatto Integration\n");
238 printf(
"Begin checking differentiation through Gauss points\n");
245 zwgj(z,w,np,alpha,beta);
246 for(n = 2; n < np-1; ++n){
247 Dgj(d,dt,z,np,alpha,beta);
249 for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
251 for(i = 0; i < np; ++i)
252 sum += fabs(
ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
255 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
256 ,alpha,beta,np,n,sum);
261 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
264 printf(
"Finished checking Gauss Jacobi differentiation\n");
267 #if GAUSS_RADAUM_DIFF
268 printf(
"Begin checking differentiation through Gauss Radau points\n");
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);
280 for(i = 0; i < np; ++i)
281 sum += fabs(
ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
284 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
285 ,alpha,beta,np,n,sum);
290 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
293 printf(
"Finished checking Gauss Radau (z=-1) differentiation\n");
296 #if GAUSS_RADAUP_DIFF
297 printf(
"Begin checking differentiation through Gauss Radau (z=1) points\n");
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);
309 for(i = 0; i < np; ++i)
310 sum += fabs(
ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
313 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
314 ,alpha,beta,np,n,sum);
319 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
322 printf(
"Finished checking Gauss Radau (z=1) differentiation\n");
325 #if GAUSS_LOBATTO_DIFF
326 printf(
"Begin checking differentiation through Gauss Lobatto points\n");
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);
338 for(i = 0; i < np; ++i)
339 sum += fabs(
ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
342 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
343 ,alpha,beta,np,n,sum);
348 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
351 printf(
"Finished checking Gauss Lobatto differentiation\n");
356 printf(
"Begin checking interpolation through Gauss points\n");
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;
369 Imgj(d,z,w,np,np,alpha,beta);
371 for(i = 0; i < np; ++i)
372 sum += fabs(
ddot(np,d+i*np,1,p,1) - pow(w[i],n));
375 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
376 ,alpha,beta,np,n,sum);
381 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
384 printf(
"Finished checking Gauss Jacobi interpolation\n");
387 #if GAUSS_RADAUM_INTERP
388 printf(
"Begin checking Interpolation through Gauss Radau (z=-1) points\n");
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;
401 Imgrjm(d,z,w,np,np,alpha,beta);
403 for(i = 0; i < np; ++i)
404 sum += fabs(
ddot(np,d+i*np,1,p,1) - pow(w[i],n));
407 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
408 ,alpha,beta,np,n,sum);
413 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
416 printf(
"Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
418 #if GAUSS_RADAUP_INTERP
419 printf(
"Begin checking Interpolation through Gauss Radau (z=1) points\n");
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;
432 Imgrjp(d,z,w,np,np,alpha,beta);
434 for(i = 0; i < np; ++i)
435 sum += fabs(
ddot(np,d+i*np,1,p,1) - pow(w[i],n));
438 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
439 ,alpha,beta,np,n,sum);
444 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
447 printf(
"Finished checking Gauss Radau (z=1) interpolation\n");
450 #if GAUSS_LOBATTO_INTERP
451 printf(
"Begin checking Interpolation through Gauss Lobatto points\n");
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;
464 Imglj(d,z,w,np,np,alpha,beta);
466 for(i = 0; i < np; ++i)
467 sum += fabs(
ddot(np,d+i*np,1,p,1) - pow(w[i],n));
470 printf(
"alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
471 ,alpha,beta,np,n,sum);
476 printf(
"finished checking all beta values for alpha = %lf\n",alpha);
479 printf(
"Finished checking Gauss Lobatto interploation\n");
483 free(z); free(w); free(p); free(d); free(dt);
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
double ddot(int, double *, int, double *, int)
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.
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.
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.
double * dvector(int, int)
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.
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
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.
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
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.
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
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.
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, .
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.