Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Polylib.cpp
Go to the documentation of this file.
1 
2 #include <stdlib.h>
3 
4 #include <stdio.h>
5 
6 #include <math.h>
7 
8 #include "Polylib.h"
9 
10 #include <float.h>
11 
12 #include <complex>
13 
15 
16 
17 /// Maximum number of iterations in polynomial defalation routine Jacobz
18 
19 #define STOP 30
20 
21 /// Precision tolerance for two points to be similar
22 
23 #define EPS 100*DBL_EPSILON
24 
25 /// return the sign(b)*a
26 
27 #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a))
28 
29 
30 
31 namespace Polylib {
32 
33  /// The following function is used to circumvent/reduce "Subtractive Cancellation"
34  /// The expression 1/dz is replaced by optinvsub(.,.)
35  /// Added on 26 April 2017
36 
37  double optdiff(double xl, double xr)
38  {
39  double m_xln, m_xrn;
40  int m_expn;
41  int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
42 
43  if (fabs(xl-xr)<1.e-4){
44 
45  m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
46  m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
47  m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
48  m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
49  m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
50 
51  return powl(10.0L,-m_digits)*(m_xln-m_xrn);
52  }else{
53  return (xl-xr);
54  }
55  }
56 
57  double laginterp(double z, int j, const double *zj, int np)
58  {
59  double temp = 1.0;
60  for (int i=0; i<np; i++)
61  {
62  if (j != i)
63  {
64  temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
65  }
66  }
67  return temp;
68  }
69 
70  /// Define whether to use polynomial deflation (1) or tridiagonal solver (0).
71 
72 #define POLYNOMIAL_DEFLATION 0
73 
74 
75 
76 #ifdef POLYNOMIAL_DEFLATION
77 
78  /// zero determination using Newton iteration with polynomial deflation
79 
80 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
81 
82 #else
83 
84  /// zero determination using eigenvalues of tridiagaonl matrix
85 
86 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
87 
88 #endif
89 
90 
91 
92 
93 
94 
95 
96  /* local functions */
97 
98  static void Jacobz (const int n, double *z, const double alpha,
99 
100  const double beta);
101 
102  // static void JacZeros (const int n, double *a, const double alpha,
103 
104  // const double beta);
105 
106  //static void TriQL (const int n, double *d, double *e);
107 
108  static void TriQL(const int, double *,double *, double **);
109 
110  double gammaF (const double);
111 
112  static void RecCoeff(const int, double *, double *,const double,
113 
114  const double);
115 
116  void JKMatrix(int, double *, double *);
117 
118  void chri1(int,double*,double*,double*,double*,double);
119 
120 
121 
122  /**
123 
124  \brief Gauss-Jacobi zeros and weights.
125 
126 
127 
128  \li Generate \a np Gauss Jacobi zeros, \a z, and weights,\a w,
129 
130  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
131 
132  \f$,
133 
134 
135 
136  \li Exact for polynomials of order \a 2np-1 or less
137 
138  */
139 
140 
141 
142  void zwgj (double *z, double *w, const int np, const double alpha,
143 
144  const double beta)
145 
146  {
147 
148  register int i;
149 
150  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
151 
152 
153 
154  jacobz (np,z,alpha,beta);
155 
156  jacobd (np,z,w,np,alpha,beta);
157 
158 
159 
160  fac = pow(two,apb + one)*gammaF(alpha + np + one)*gammaF(beta + np + one);
161 
162  fac /= gammaF(np + one)*gammaF(apb + np + one);
163 
164 
165 
166  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
167 
168 
169 
170  return;
171 
172  }
173 
174 
175 
176 
177 
178  /**
179 
180  \brief Gauss-Radau-Jacobi zeros and weights with end point at \a z=-1.
181 
182 
183 
184  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
185 
186  associated with the polynomial \f$(1+z) P^{\alpha,\beta+1}_{np-1}(z)
187 
188  \f$.
189 
190 
191 
192  \li Exact for polynomials of order \a 2np-2 or less
193 
194  */
195 
196 
197 
198  void zwgrjm(double *z, double *w, const int np, const double alpha,
199 
200  const double beta)
201 
202  {
203 
204 
205 
206  if(np == 1){
207 
208  z[0] = 0.0;
209 
210  w[0] = 2.0;
211 
212  }
213 
214  else{
215 
216  register int i;
217 
218  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
219 
220 
221 
222  z[0] = -one;
223 
224  jacobz (np-1,z+1,alpha,beta+1);
225 
226  jacobfd (np,z,w,NULL,np-1,alpha,beta);
227 
228 
229 
230  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
231 
232  fac /= gammaF(np)*(beta + np)*gammaF(apb + np + 1);
233 
234 
235 
236  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
237 
238  w[0] *= (beta + one);
239 
240  }
241 
242 
243 
244  return;
245 
246  }
247 
248 
249 
250 
251 
252  /**
253 
254  \brief Gauss-Radau-Jacobi zeros and weights with end point at \a z=1
255 
256 
257 
258 
259 
260  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
261 
262  associated with the polynomial \f$(1-z) P^{\alpha+1,\beta}_{np-1}(z)
263 
264  \f$.
265 
266 
267 
268  \li Exact for polynomials of order \a 2np-2 or less
269 
270  */
271 
272 
273 
274  void zwgrjp(double *z, double *w, const int np, const double alpha,
275 
276  const double beta)
277 
278  {
279 
280 
281 
282  if(np == 1){
283 
284  z[0] = 0.0;
285 
286  w[0] = 2.0;
287 
288  }
289 
290  else{
291 
292  register int i;
293 
294  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
295 
296 
297 
298  jacobz (np-1,z,alpha+1,beta);
299 
300  z[np-1] = one;
301 
302  jacobfd (np,z,w,NULL,np-1,alpha,beta);
303 
304 
305 
306  fac = pow(two,apb)*gammaF(alpha + np)*gammaF(beta + np);
307 
308  fac /= gammaF(np)*(alpha + np)*gammaF(apb + np + 1);
309 
310 
311 
312  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
313 
314  w[np-1] *= (alpha + one);
315 
316  }
317 
318 
319 
320  return;
321 
322  }
323 
324 
325 
326 
327 
328  /**
329 
330  \brief Gauss-Lobatto-Jacobi zeros and weights with end point at \a z=-1,\a 1
331 
332 
333 
334 
335 
336  \li Generate \a np Gauss-Lobatto-Jacobi points, \a z, and weights, \a w,
337 
338  associated with polynomial \f$ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \f$
339 
340  \li Exact for polynomials of order \a 2np-3 or less
341 
342  */
343 
344 
345 
346  void zwglj(double *z, double *w, const int np, const double alpha,
347 
348  const double beta)
349 
350  {
351 
352 
353 
354  if( np == 1 ){
355 
356  z[0] = 0.0;
357 
358  w[0] = 2.0;
359 
360  }
361 
362  else{
363 
364  register int i;
365 
366  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
367 
368 
369 
370  z[0] = -one;
371 
372  z[np-1] = one;
373 
374  jacobz (np-2,z + 1,alpha + one,beta + one);
375 
376  jacobfd (np,z,w,NULL,np-1,alpha,beta);
377 
378 
379 
380  fac = pow(two,apb + 1)*gammaF(alpha + np)*gammaF(beta + np);
381 
382  fac /= (np-1)*gammaF(np)*gammaF(alpha + beta + np + one);
383 
384 
385 
386  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
387 
388  w[0] *= (beta + one);
389 
390  w[np-1] *= (alpha + one);
391 
392  }
393 
394 
395 
396  return;
397 
398  }
399 
400 
401 
402  /**
403 
404  \brief Gauss-Kronrod-Jacobi zeros and weights.
405 
406 
407 
408  \li Generate \a npt=2*np+1 Gauss-Kronrod Jacobi zeros, \a z, and weights,\a w,
409 
410  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
411 
412  \f$,
413 
414 
415 
416  \li Exact for polynomials of order \a 3np+1 or less
417 
418  */
419 
420  void zwgk(double *z, double *w, const int npt , const double alpha,
421 
422  const double beta)
423 
424  {
425 
426 
427 
428  int np = (npt-1)/2;
429 
430 
431 
432  int i,j;
433 
434 
435 
436  // number of kronrod points associated with the np gauss rule
437 
438  int kpoints = 2*np + 1;
439 
440 
441 
442  // Define the number of required recurrence coefficents
443 
444  int ncoeffs = (int)floor(3.0*(np+1)/2);
445 
446 
447 
448  // Define arrays for the recurrence coefficients
449 
450  // We will use these arrays for the Kronrod results too, hence the
451 
452  // reason for the size of the arrays
453 
454  double *a = new double[kpoints];
455 
456  double *b = new double[kpoints];
457 
458 
459 
460  // Initialize a and b to zero
461 
462  for(i = 0; i < kpoints; i++)
463 
464  {
465 
466  a[i] = 0.0;
467 
468  b[i] = 0.0;
469 
470  }
471 
472 
473 
474  // Call the routine to calculate the recurrence coefficients
475 
476  RecCoeff(ncoeffs,a,b,alpha,beta);
477 
478 
479 
480  // Call the routine to caluclate the jacobi-Kronrod matrix
481 
482  JKMatrix(np,a,b);
483 
484 
485 
486  // Set up the identity matrix
487 
488  double** zmatrix = new double*[kpoints];
489 
490  for(i = 0; i < kpoints; i++)
491 
492  {
493 
494  zmatrix[i] = new double[kpoints];
495 
496  for(j = 0; j < kpoints; j++)
497 
498  {
499 
500  zmatrix[i][j] = 0.0;
501 
502  }
503 
504  }
505 
506  for(i = 0; i < kpoints; i++)
507 
508  {
509 
510  zmatrix[i][i] = 1.0;
511 
512  }
513 
514 
515 
516  // Calculte the points and weights
517 
518  TriQL(kpoints, a, b, zmatrix);
519 
520 
521 
522  for(i = 0; i < kpoints; i++)
523 
524  {
525 
526  z[i] = a[i];
527 
528  w[i] = b[i];
529 
530  }
531 
532 
533 
534 
535 
536  }
537 
538 
539 
540  /**
541 
542  \brief Gauss-Radau-Kronrod-Jacobi zeros and weights.
543 
544 
545 
546  \li Generate \a npt=2*np Radau-Kronrod Jacobi zeros, \a z, and weights,\a w,
547 
548  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
549 
550  \f$,
551 
552  */
553 
554  void zwrk(double* z, double* w, const int npt ,const double alpha,
555 
556  const double beta)
557 
558  {
559 
560 
561 
562  int np = npt/2;
563 
564 
565 
566  if(np < 2)
567 
568  {
569 
570  fprintf(stderr,"too few points in formula\n");
571 
572  return;
573 
574  }
575 
576 
577 
578  double end0 = -1;
579 
580 
581 
582  int i,j;
583 
584 
585 
586  // number of kronrod points associated with the np gauss rule
587 
588  int kpoints = 2*np;
589 
590 
591 
592  // Define the number of required recurrence coefficents
593 
594  int ncoeffs = (int)ceil(3.0*np/2);
595 
596 
597 
598  // Define arrays for the recurrence coefficients
599 
600  double *a = new double[ncoeffs+1];
601 
602  double *b = new double[ncoeffs+1];
603 
604 
605 
606  // Initialize a and b to zero
607 
608  for(i = 0; i < ncoeffs+1; i++)
609 
610  {
611 
612  a[i] = 0.0;
613 
614  b[i] = 0.0;
615 
616  }
617 
618 
619 
620  // Call the routine to calculate the recurrence coefficients
621 
622  RecCoeff(ncoeffs,a,b,alpha,beta);
623 
624 
625 
626  double* a0 = new double[ncoeffs];
627 
628  double* b0 = new double[ncoeffs];
629 
630 
631 
632  chri1(ncoeffs,a,b,a0,b0,end0);
633 
634 
635 
636  double s = b0[0]/fabs(b0[0]);
637 
638  b0[0] = s*b0[0];
639 
640 
641 
642  // Finding the 2*np-1 gauss-kronrod points
643 
644  double* z1 = new double[2*np-1];
645 
646  double* w1 = new double[2*np-1];
647 
648  for(i = 0; i < ncoeffs; i++)
649 
650  {
651 
652  z1[i] = a0[i];
653 
654  w1[i] = b0[i];
655 
656  }
657 
658  JKMatrix(np-1,z1,w1);
659 
660  // Set up the identity matrix
661 
662  double** zmatrix = new double*[2*np-1];
663 
664  for(i = 0; i < 2*np-1; i++)
665 
666  {
667 
668  zmatrix[i] = new double[2*np-1];
669 
670  for(j = 0; j < 2*np-1; j++)
671 
672  {
673 
674  zmatrix[i][j] = 0.0;
675 
676  }
677 
678  }
679 
680  for(i = 0; i < 2*np-1; i++)
681 
682  {
683 
684  zmatrix[i][i] = 1.0;
685 
686  }
687 
688 
689 
690  // Calculate the points and weights
691 
692  TriQL(2*np-1, z1, w1, zmatrix);
693 
694 
695 
696  double sumW1 = 0.0;
697 
698  for(i = 0; i < 2*np-1; i++)
699 
700  {
701 
702  w1[i] = s*w1[i]/(z1[i]-end0);
703 
704  sumW1 += w1[i];
705 
706  }
707 
708 
709 
710  z[0] = end0;
711 
712  w[0] = b[0]- sumW1;
713 
714  for(i = 1; i < kpoints; i++)
715 
716  {
717 
718  z[i] = z1[i-1];
719 
720  w[i] = w1[i-1];
721 
722  }
723 
724 
725 
726 
727 
728  }
729 
730 
731 
732  /**
733 
734  \brief Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
735 
736 
737 
738  \li Generate \a npt=2*np-1 Lobatto-Kronrod Jacobi zeros, \a z, and weights,\a w,
739 
740  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
741 
742  \f$,
743 
744  */
745 
746  void zwlk(double* z, double* w, const int npt,
747 
748  const double alpha, const double beta)
749 
750  {
751 
752 
753 
754  int np = (npt+1)/2;
755 
756 
757 
758  if(np < 4)
759 
760  {
761 
762  fprintf (stderr,"too few points in formula\n");
763 
764  return;
765 
766  }
767 
768 
769 
770  double endl = -1;
771 
772  double endr = 1;
773 
774  int i,j;
775 
776 
777 
778  // number of kronrod points associated with the np gauss rule
779 
780  int kpoints = 2*np-1;
781 
782 
783 
784  // Define the number of required recurrence coefficents
785 
786  int ncoeffs = (int)ceil(3.0*np/2)-1;
787 
788 
789 
790  // Define arrays for the recurrence coefficients
791 
792  double *a = new double[ncoeffs+1];
793 
794  double *b = new double[ncoeffs+1];
795 
796 
797 
798  // Initialize a and b to zero
799 
800  for(i = 0; i < ncoeffs+1; i++)
801 
802  {
803 
804  a[i] = 0.0;
805 
806  b[i] = 0.0;
807 
808  }
809 
810 
811 
812  // Call the routine to calculate the recurrence coefficients
813 
814  RecCoeff(ncoeffs,a,b,alpha,beta);
815 
816 
817 
818 
819 
820  double* a0 = new double[ncoeffs];
821 
822  double* b0 = new double[ncoeffs];
823 
824 
825 
826  chri1(ncoeffs,a,b,a0,b0,endl);
827 
828 
829 
830  double* a1 = new double[ncoeffs-1];
831 
832  double* b1 = new double[ncoeffs-1];
833 
834 
835 
836  chri1(ncoeffs-1,a0,b0,a1,b1,endr);
837 
838 
839 
840 
841 
842  double s = b1[0]/fabs(b1[0]);
843 
844  b1[0] = s*b1[0];
845 
846 
847 
848  // Finding the 2*np-1 gauss-kronrod points
849 
850  double* z1 = new double[2*np-3];
851 
852  double* w1 = new double[2*np-3];
853 
854  for(i = 0; i < ncoeffs; i++)
855 
856  {
857 
858  z1[i] = a1[i];
859 
860  w1[i] = b1[i];
861 
862  }
863 
864  JKMatrix(np-2,z1,w1);
865 
866  // Set up the identity matrix
867 
868  double** zmatrix = new double*[2*np-3];
869 
870  for(i = 0; i < 2*np-3; i++)
871 
872  {
873 
874  zmatrix[i] = new double[2*np-3];
875 
876  for(j = 0; j < 2*np-3; j++)
877 
878  {
879 
880  zmatrix[i][j] = 0.0;
881 
882  }
883 
884  }
885 
886  for(i = 0; i < 2*np-3; i++)
887 
888  {
889 
890  zmatrix[i][i] = 1.0;
891 
892  }
893 
894 
895 
896  // Calculate the points and weights
897 
898  TriQL(2*np-3, z1, w1, zmatrix);
899 
900 
901 
902  double sumW1 = 0.0;
903 
904  double sumW1Z1 = 0.0;
905 
906  for(i = 0; i < 2*np-3; i++)
907 
908  {
909 
910  w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
911 
912  sumW1 += w1[i];
913 
914  sumW1Z1 += z1[i]*w1[i];
915 
916  }
917 
918 
919 
920  double c0 = b[0]-sumW1;
921 
922  double c1 = a[0]*b[0]-sumW1Z1;
923 
924 
925 
926  z[0] = endl;
927 
928  z[2*np-2] = endr;
929 
930  w[0] = (c0*endr-c1)/(endr-endl);
931 
932  w[2*np-2] = (c1-c0*endl)/(endr-endl);
933 
934 
935 
936  for(i = 1; i < kpoints-1; i++)
937 
938  {
939 
940  z[i] = z1[i-1];
941 
942  w[i] = w1[i-1];
943 
944  }
945 
946  }
947 
948 
949 
950  /**
951 
952  \brief Compute the Derivative Matrix and its transpose associated
953 
954  with the Gauss-Jacobi zeros.
955 
956 
957 
958  \li Compute the derivative matrix, \a d, and its transpose, \a dt,
959 
960  associated with the n_th order Lagrangian interpolants through the
961 
962  \a np Gauss-Jacobi points \a z such that \n
963 
964  \f$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
965 
966 
967 
968  */
969 
970 
971 
972  void Dgj(double *D, const double *z, const int np, const double alpha,
973 
974  const double beta)
975 
976  {
977 
978 
979 
980  double one = 1.0, two = 2.0;
981 
982 
983 
984  if (np <= 0){
985 
986  D[0] = 0.0;
987 
988  }
989 
990  else{
991 
992  register int i,j;
993 
994  double *pd;
995 
996 
997 
998  pd = (double *)malloc(np*sizeof(double));
999 
1000  jacobd(np,z,pd,np,alpha,beta);
1001 
1002 
1003 
1004  for (i = 0; i < np; i++){
1005 
1006  for (j = 0; j < np; j++){
1007 
1008 
1009 
1010  if (i != j)
1011 
1012  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1013 
1014  else
1015 
1016  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
1017 
1018  (two*(one - z[j]*z[j]));
1019 
1020  }
1021 
1022  }
1023 
1024  free(pd);
1025 
1026  }
1027 
1028  return;
1029 
1030  }
1031 
1032 
1033 
1034 
1035 
1036  /**
1037 
1038  \brief Compute the Derivative Matrix and its transpose associated
1039 
1040  with the Gauss-Radau-Jacobi zeros with a zero at \a z=-1.
1041 
1042 
1043 
1044  \li Compute the derivative matrix, \a d, associated with the n_th
1045 
1046  order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
1047 
1048  points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
1049 
1050  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
1051 
1052 
1053 
1054  */
1055 
1056 
1057 
1058  void Dgrjm(double *D, const double *z, const int np, const double alpha,
1059 
1060  const double beta)
1061 
1062  {
1063 
1064 
1065 
1066  if (np <= 0){
1067 
1068  D[0] = 0.0;
1069 
1070  }
1071 
1072  else{
1073 
1074  register int i, j;
1075 
1076  double one = 1.0, two = 2.0;
1077 
1078  double *pd;
1079 
1080 
1081 
1082  pd = (double *)malloc(np*sizeof(double));
1083 
1084 
1085 
1086  pd[0] = pow(-one,np-1)*gammaF(np+beta+one);
1087 
1088  pd[0] /= gammaF(np)*gammaF(beta+two);
1089 
1090  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
1091 
1092  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
1093 
1094 
1095 
1096  for (i = 0; i < np; i++) {
1097 
1098  for (j = 0; j < np; j++){
1099 
1100  if (i != j)
1101 
1102  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1103 
1104  else {
1105 
1106  if(j == 0)
1107 
1108  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
1109 
1110  (two*(beta + two));
1111 
1112  else
1113 
1114  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
1115 
1116  (two*(one - z[j]*z[j]));
1117 
1118  }
1119 
1120  }
1121 
1122  }
1123 
1124  free(pd);
1125 
1126  }
1127 
1128 
1129 
1130  return;
1131 
1132  }
1133 
1134 
1135 
1136 
1137 
1138  /**
1139 
1140  \brief Compute the Derivative Matrix associated with the
1141 
1142  Gauss-Radau-Jacobi zeros with a zero at \a z=1.
1143 
1144 
1145 
1146  \li Compute the derivative matrix, \a d, associated with the n_th
1147 
1148  order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
1149 
1150  points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
1151 
1152  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
1153 
1154  */
1155 
1156 
1157 
1158  void Dgrjp(double *D, const double *z, const int np, const double alpha,
1159 
1160  const double beta)
1161 
1162  {
1163 
1164 
1165 
1166  if (np <= 0){
1167 
1168  D[0] = 0.0;
1169 
1170  }
1171 
1172  else{
1173 
1174  register int i, j;
1175 
1176  double one = 1.0, two = 2.0;
1177 
1178  double *pd;
1179 
1180 
1181 
1182  pd = (double *)malloc(np*sizeof(double));
1183 
1184 
1185 
1186 
1187 
1188  jacobd(np-1,z,pd,np-1,alpha+1,beta);
1189 
1190  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
1191 
1192  pd[np-1] = -gammaF(np+alpha+one);
1193 
1194  pd[np-1] /= gammaF(np)*gammaF(alpha+two);
1195 
1196 
1197 
1198  for (i = 0; i < np; i++) {
1199 
1200  for (j = 0; j < np; j++){
1201 
1202  if (i != j)
1203 
1204  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1205 
1206  else {
1207 
1208  if(j == np-1)
1209 
1210  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
1211 
1212  (two*(alpha + two));
1213 
1214  else
1215 
1216  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
1217 
1218  (two*(one - z[j]*z[j]));
1219 
1220  }
1221 
1222  }
1223 
1224  }
1225 
1226  free(pd);
1227 
1228  }
1229 
1230 
1231 
1232  return;
1233 
1234  }
1235 
1236 
1237 
1238  /**
1239 
1240  \brief Compute the Derivative Matrix associated with the
1241 
1242  Gauss-Lobatto-Jacobi zeros.
1243 
1244 
1245 
1246  \li Compute the derivative matrix, \a d, associated with the n_th
1247 
1248  order Lagrange interpolants through the \a np
1249 
1250  Gauss-Lobatto-Jacobi points \a z such that \n \f$
1251 
1252  \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
1253 
1254 
1255 
1256  */
1257 
1258 
1259 
1260  void Dglj(double *D, const double *z, const int np, const double alpha,
1261 
1262  const double beta)
1263 
1264  {
1265 
1266 
1267 
1268  if (np <= 0){
1269 
1270  D[0] = 0.0;
1271 
1272  }
1273 
1274  else{
1275 
1276  register int i, j;
1277 
1278  double one = 1.0, two = 2.0;
1279 
1280  double *pd;
1281 
1282 
1283 
1284  pd = (double *)malloc(np*sizeof(double));
1285 
1286 
1287 
1288  pd[0] = two*pow(-one,np)*gammaF(np + beta);
1289 
1290  pd[0] /= gammaF(np - one)*gammaF(beta + two);
1291 
1292  jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
1293 
1294  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
1295 
1296  pd[np-1] = -two*gammaF(np + alpha);
1297 
1298  pd[np-1] /= gammaF(np - one)*gammaF(alpha + two);
1299 
1300 
1301 
1302  for (i = 0; i < np; i++) {
1303 
1304  for (j = 0; j < np; j++){
1305 
1306  if (i != j)
1307 
1308  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1309 
1310  else {
1311 
1312  if (j == 0)
1313 
1314  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
1315 
1316  else if (j == np-1)
1317 
1318  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
1319 
1320  else
1321 
1322  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
1323 
1324  (two*(one - z[j]*z[j]));
1325 
1326  }
1327 
1328  }
1329 
1330  }
1331 
1332  free(pd);
1333 
1334  }
1335 
1336 
1337 
1338  return;
1339 
1340  }
1341 
1342 
1343 
1344 
1345 
1346  /**
1347 
1348  \brief Compute the value of the \a i th Lagrangian interpolant through
1349 
1350  the \a np Gauss-Jacobi points \a zgj at the arbitrary location \a z.
1351 
1352 
1353 
1354  \li \f$ -1 \leq z \leq 1 \f$
1355 
1356 
1357 
1358  \li Uses the defintion of the Lagrangian interpolant:\n
1359 
1360  %
1361 
1362  \f$ \begin{array}{rcl}
1363 
1364  h_j(z) = \left\{ \begin{array}{ll}
1365 
1366  \displaystyle \frac{P_{np}^{\alpha,\beta}(z)}
1367 
1368  {[P_{np}^{\alpha,\beta}(z_j)]^\prime
1369 
1370  (z-z_j)} & \mbox{if $z \ne z_j$}\\
1371 
1372  & \\
1373 
1374  1 & \mbox{if $z=z_j$}
1375 
1376  \end{array}
1377 
1378  \right.
1379 
1380  \end{array} \f$
1381 
1382  */
1383 
1384 
1385 
1386  double hgj (const int i, const double z, const double *zgj,
1387 
1388  const int np, const double alpha, const double beta)
1389 
1390  {
1391 
1392  double zi, dz;
1393 
1394  zi = *(zgj+i);
1395 
1396  dz = z-zi;
1397 
1398  if (fabs(dz) < EPS) return 1.0;
1399 
1400  return laginterp(z, i, zgj, np);
1401 
1402  }
1403 
1404 
1405 
1406  /**
1407 
1408  \brief Compute the value of the \a i th Lagrangian interpolant through the
1409 
1410  \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
1411 
1412  \a z. This routine assumes \a zgrj includes the point \a -1.
1413 
1414 
1415 
1416  \li \f$ -1 \leq z \leq 1 \f$
1417 
1418 
1419 
1420  \li Uses the defintion of the Lagrangian interpolant:\n
1421 
1422  %
1423 
1424  \f$ \begin{array}{rcl}
1425 
1426  h_j(z) = \left\{ \begin{array}{ll}
1427 
1428  \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)}
1429 
1430  {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime +
1431 
1432  P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
1433 
1434  & \\
1435 
1436  1 & \mbox{if $z=z_j$}
1437 
1438  \end{array}
1439 
1440  \right.
1441 
1442  \end{array} \f$
1443 
1444  */
1445 
1446 
1447 
1448  double hgrjm (const int i, const double z, const double *zgrj, const int np,
1449 
1450  const double alpha, const double beta)
1451 
1452  {
1453 
1454 
1455  double zi, dz;
1456 
1457  zi = *(zgrj+i);
1458 
1459  dz = z-zi;
1460 
1461  if (fabs(dz) < EPS) return 1.0;
1462 
1463  return laginterp(z, i, zgrj, np);
1464 
1465  }
1466 
1467 
1468 
1469 
1470 
1471  /**
1472 
1473  \brief Compute the value of the \a i th Lagrangian interpolant through the
1474 
1475  \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
1476 
1477  \a z. This routine assumes \a zgrj includes the point \a +1.
1478 
1479 
1480 
1481  \li \f$ -1 \leq z \leq 1 \f$
1482 
1483 
1484 
1485  \li Uses the defintion of the Lagrangian interpolant:\n
1486 
1487  %
1488 
1489  \f$ \begin{array}{rcl}
1490 
1491  h_j(z) = \left\{ \begin{array}{ll}
1492 
1493  \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)}
1494 
1495  {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime -
1496 
1497  P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
1498 
1499  & \\
1500 
1501  1 & \mbox{if $z=z_j$}
1502 
1503  \end{array}
1504 
1505  \right.
1506 
1507  \end{array} \f$
1508 
1509  */
1510 
1511 
1512 
1513  double hgrjp (const int i, const double z, const double *zgrj, const int np,
1514 
1515  const double alpha, const double beta)
1516 
1517  {
1518 
1519 
1520  double zi, dz;
1521 
1522  zi = *(zgrj+i);
1523 
1524  dz = z-zi;
1525 
1526  if (fabs(dz) < EPS) return 1.0;
1527 
1528  return laginterp(z, i, zgrj, np);
1529 
1530  }
1531 
1532 
1533 
1534 
1535 
1536  /**
1537 
1538  \brief Compute the value of the \a i th Lagrangian interpolant through the
1539 
1540  \a np Gauss-Lobatto-Jacobi points \a zgrj at the arbitrary location
1541 
1542  \a z.
1543 
1544 
1545 
1546  \li \f$ -1 \leq z \leq 1 \f$
1547 
1548 
1549 
1550  \li Uses the defintion of the Lagrangian interpolant:\n
1551 
1552  %
1553 
1554  \f$ \begin{array}{rcl}
1555 
1556  h_j(z) = \left\{ \begin{array}{ll}
1557 
1558  \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)}
1559 
1560  {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime -
1561 
1562  2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\
1563 
1564  & \\
1565 
1566  1 & \mbox{if $z=z_j$}
1567 
1568  \end{array}
1569 
1570  \right.
1571 
1572  \end{array} \f$
1573 
1574  */
1575 
1576 
1577 
1578  double hglj (const int i, const double z, const double *zglj, const int np,
1579 
1580  const double alpha, const double beta)
1581 
1582  {
1583 
1584  double zi, dz;
1585 
1586  zi = *(zglj+i);
1587 
1588  dz = z-zi;
1589 
1590  if (fabs(dz) < EPS) return 1.0;
1591 
1592  return laginterp(z, i, zglj, np);
1593 
1594  }
1595 
1596 
1597 
1598 
1599 
1600  /**
1601 
1602  \brief Interpolation Operator from Gauss-Jacobi points to an
1603 
1604  arbitrary distribution at points \a zm
1605 
1606 
1607 
1608  \li Computes the one-dimensional interpolation matrix, \a im, to
1609 
1610  interpolate a function from at Gauss-Jacobi distribution of \a nz
1611 
1612  zeros \a zgrj to an arbitrary distribution of \a mz points \a zm, i.e.\n
1613 
1614  \f$
1615 
1616  u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j])
1617 
1618  \f$
1619 
1620 
1621 
1622  */
1623 
1624 
1625 
1626  void Imgj(double *im, const double *zgj, const double *zm, const int nz,
1627 
1628  const int mz,const double alpha, const double beta){
1629 
1630  double zp;
1631 
1632  register int i, j;
1633 
1634 
1635 
1636  for (i = 0; i < nz; ++i) {
1637 
1638  for (j = 0; j < mz; ++j)
1639 
1640  {
1641 
1642  zp = zm[j];
1643 
1644  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
1645 
1646  }
1647 
1648  }
1649 
1650 
1651 
1652  return;
1653 
1654  }
1655 
1656 
1657 
1658  /**
1659 
1660  \brief Interpolation Operator from Gauss-Radau-Jacobi points
1661 
1662  (including \a z=-1) to an arbitrary distrubtion at points \a zm
1663 
1664 
1665 
1666  \li Computes the one-dimensional interpolation matrix, \a im, to
1667 
1668  interpolate a function from at Gauss-Radau-Jacobi distribution of
1669 
1670  \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
1671 
1672  distribution of \a mz points \a zm, i.e.
1673 
1674  \n
1675 
1676  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
1677 
1678 
1679 
1680  */
1681 
1682 
1683 
1684  void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz,
1685 
1686  const int mz, const double alpha, const double beta)
1687 
1688  {
1689 
1690  double zp;
1691 
1692  register int i, j;
1693 
1694 
1695 
1696  for (i = 0; i < nz; i++) {
1697 
1698  for (j = 0; j < mz; j++)
1699 
1700  {
1701 
1702  zp = zm[j];
1703 
1704  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1705 
1706  }
1707 
1708  }
1709 
1710 
1711 
1712  return;
1713 
1714  }
1715 
1716 
1717 
1718  /**
1719 
1720  \brief Interpolation Operator from Gauss-Radau-Jacobi points
1721 
1722  (including \a z=1) to an arbitrary distrubtion at points \a zm
1723 
1724 
1725 
1726  \li Computes the one-dimensional interpolation matrix, \a im, to
1727 
1728  interpolate a function from at Gauss-Radau-Jacobi distribution of
1729 
1730  \a nz zeros \a zgrj (where \a zgrj[nz-1]=1) to an arbitrary
1731 
1732  distribution of \a mz points \a zm, i.e.
1733 
1734  \n
1735 
1736  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
1737 
1738 
1739 
1740  */
1741 
1742 
1743 
1744  void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz,
1745 
1746  const int mz,const double alpha, const double beta)
1747 
1748  {
1749 
1750  double zp;
1751 
1752  register int i, j;
1753 
1754 
1755 
1756  for (i = 0; i < nz; i++) {
1757 
1758  for (j = 0; j < mz; j++)
1759 
1760  {
1761 
1762  zp = zm[j];
1763 
1764  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1765 
1766  }
1767 
1768  }
1769 
1770 
1771 
1772  return;
1773 
1774  }
1775 
1776 
1777 
1778 
1779 
1780  /**
1781 
1782  \brief Interpolation Operator from Gauss-Lobatto-Jacobi points
1783 
1784  to an arbitrary distrubtion at points \a zm
1785 
1786 
1787 
1788  \li Computes the one-dimensional interpolation matrix, \a im, to
1789 
1790  interpolate a function from at Gauss-Lobatto-Jacobi distribution of
1791 
1792  \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
1793 
1794  distribution of \a mz points \a zm, i.e.
1795 
1796  \n
1797 
1798  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
1799 
1800 
1801 
1802  */
1803 
1804 
1805 
1806  void Imglj(double *im, const double *zglj, const double *zm, const int nz,
1807 
1808  const int mz, const double alpha, const double beta)
1809 
1810  {
1811 
1812  double zp;
1813 
1814  register int i, j;
1815 
1816 
1817 
1818  for (i = 0; i < nz; i++) {
1819 
1820  for (j = 0; j < mz; j++)
1821 
1822  {
1823 
1824  zp = zm[j];
1825 
1826  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
1827 
1828  }
1829 
1830  }
1831 
1832 
1833 
1834  return;
1835 
1836  }
1837 
1838 
1839 
1840  /**
1841 
1842  \brief Routine to calculate Jacobi polynomials, \f$
1843 
1844  P^{\alpha,\beta}_n(z) \f$, and their first derivative, \f$
1845 
1846  \frac{d}{dz} P^{\alpha,\beta}_n(z) \f$.
1847 
1848 
1849 
1850  \li This function returns the vectors \a poly_in and \a poly_d
1851 
1852  containing the value of the \f$ n^th \f$ order Jacobi polynomial
1853 
1854  \f$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \f$ and its
1855 
1856  derivative at the \a np points in \a z[i]
1857 
1858 
1859 
1860  - If \a poly_in = NULL then only calculate derivatice
1861 
1862 
1863 
1864  - If \a polyd = NULL then only calculate polynomial
1865 
1866 
1867 
1868  - To calculate the polynomial this routine uses the recursion
1869 
1870  relationship (see appendix A ref [4]) :
1871 
1872  \f$ \begin{array}{rcl}
1873 
1874  P^{\alpha,\beta}_0(z) &=& 1 \\
1875 
1876  P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\
1877 
1878  a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z)
1879 
1880  P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\
1881 
1882  a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\
1883 
1884  a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\
1885 
1886  a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1)
1887 
1888  (2n + \alpha + \beta + 2) \\
1889 
1890  a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2)
1891 
1892  \end{array} \f$
1893 
1894 
1895 
1896  - To calculate the derivative of the polynomial this routine uses
1897 
1898  the relationship (see appendix A ref [4]) :
1899 
1900  \f$ \begin{array}{rcl}
1901 
1902  b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z)
1903 
1904  + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\
1905 
1906  b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\
1907 
1908  b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\
1909 
1910  b^3_n(z) &=& 2(n+\alpha)(n+\beta)
1911 
1912  \end{array} \f$
1913 
1914 
1915 
1916  - Note the derivative from this routine is only valid for -1 < \a z < 1.
1917 
1918  */
1919 
1920  void jacobfd(const int np, const double *z, double *poly_in, double *polyd,
1921 
1922  const int n, const double alpha, const double beta){
1923 
1924  register int i;
1925 
1926  double zero = 0.0, one = 1.0, two = 2.0;
1927 
1928 
1929 
1930  if(!np)
1931 
1932  return;
1933 
1934 
1935 
1936  if(n == 0){
1937 
1938  if(poly_in)
1939 
1940  for(i = 0; i < np; ++i)
1941 
1942  poly_in[i] = one;
1943 
1944  if(polyd)
1945 
1946  for(i = 0; i < np; ++i)
1947 
1948  polyd[i] = zero;
1949 
1950  }
1951 
1952  else if (n == 1){
1953 
1954  if(poly_in)
1955 
1956  for(i = 0; i < np; ++i)
1957 
1958  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1959 
1960  if(polyd)
1961 
1962  for(i = 0; i < np; ++i)
1963 
1964  polyd[i] = 0.5*(alpha + beta + two);
1965 
1966  }
1967 
1968  else{
1969 
1970  register int k;
1971 
1972  double a1,a2,a3,a4;
1973 
1974  double two = 2.0, apb = alpha + beta;
1975 
1976  double *poly, *polyn1,*polyn2;
1977 
1978 
1979 
1980  if(poly_in){ // switch for case of no poynomial function return
1981 
1982  polyn1 = (double *)malloc(2*np*sizeof(double));
1983 
1984  polyn2 = polyn1+np;
1985 
1986  poly = poly_in;
1987 
1988  }
1989 
1990  else{
1991 
1992  polyn1 = (double *)malloc(3*np*sizeof(double));
1993 
1994  polyn2 = polyn1+np;
1995 
1996  poly = polyn2+np;
1997 
1998  }
1999 
2000 
2001 
2002  for(i = 0; i < np; ++i){
2003 
2004  polyn2[i] = one;
2005 
2006  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
2007 
2008  }
2009 
2010 
2011 
2012  for(k = 2; k <= n; ++k){
2013 
2014  a1 = two*k*(k + apb)*(two*k + apb - two);
2015 
2016  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
2017 
2018  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
2019 
2020  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
2021 
2022 
2023 
2024  a2 /= a1;
2025 
2026  a3 /= a1;
2027 
2028  a4 /= a1;
2029 
2030 
2031 
2032  for(i = 0; i < np; ++i){
2033 
2034  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
2035 
2036  polyn2[i] = polyn1[i];
2037 
2038  polyn1[i] = poly [i];
2039 
2040  }
2041 
2042  }
2043 
2044 
2045 
2046  if(polyd){
2047 
2048  a1 = n*(alpha - beta);
2049 
2050  a2 = n*(two*n + alpha + beta);
2051 
2052  a3 = two*(n + alpha)*(n + beta);
2053 
2054  a4 = (two*n + alpha + beta);
2055 
2056  a1 /= a4; a2 /= a4; a3 /= a4;
2057 
2058 
2059 
2060  // note polyn2 points to polyn1 at end of poly iterations
2061 
2062  for(i = 0; i < np; ++i){
2063 
2064  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
2065 
2066  polyd[i] /= (one - z[i]*z[i]);
2067 
2068  }
2069 
2070  }
2071 
2072 
2073 
2074  free(polyn1);
2075 
2076  }
2077 
2078 
2079 
2080  return;
2081 
2082  }
2083 
2084 
2085 
2086 
2087 
2088  /**
2089 
2090  \brief Calculate the derivative of Jacobi polynomials
2091 
2092 
2093 
2094  \li Generates a vector \a poly of values of the derivative of the
2095 
2096  \a n th order Jacobi polynomial \f$ P^(\alpha,\beta)_n(z)\f$ at the
2097 
2098  \a np points \a z.
2099 
2100 
2101 
2102  \li To do this we have used the relation
2103 
2104  \n
2105 
2106  \f$ \frac{d}{dz} P^{\alpha,\beta}_n(z)
2107 
2108  = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) \f$
2109 
2110 
2111 
2112  \li This formulation is valid for \f$ -1 \leq z \leq 1 \f$
2113 
2114 
2115 
2116  */
2117 
2118 
2119 
2120  void jacobd(const int np, const double *z, double *polyd, const int n,
2121 
2122  const double alpha, const double beta)
2123 
2124  {
2125 
2126  register int i;
2127 
2128  double one = 1.0;
2129 
2130  if(n == 0)
2131 
2132  for(i = 0; i < np; ++i) polyd[i] = 0.0;
2133 
2134  else{
2135 
2136  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
2137 
2138  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2139 
2140  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
2141 
2142  }
2143 
2144  return;
2145 
2146  }
2147 
2148 
2149 
2150 
2151 
2152  /**
2153 
2154  \brief Calculate the Gamma function , \f$ \Gamma(n)\f$, for integer
2155 
2156  values and halves.
2157 
2158 
2159 
2160  Determine the value of \f$\Gamma(n)\f$ using:
2161 
2162 
2163 
2164  \f$ \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$
2165 
2166 
2167 
2168  where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
2169 
2170  */
2171 
2172 
2173 
2174  double gammaF(const double x){
2175 
2176  double gamma = 1.0;
2177 
2178 
2179 
2180  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2181 
2182  else if (!x) return gamma;
2183 
2184  else if ((x-(int)x) == 0.5){
2185 
2186  int n = (int) x;
2187 
2188  double tmp = x;
2189 
2190 
2191 
2192  gamma = sqrt(M_PI);
2193 
2194  while(n--){
2195 
2196  tmp -= 1.0;
2197 
2198  gamma *= tmp;
2199 
2200  }
2201 
2202  }
2203 
2204  else if ((x-(int)x) == 0.0){
2205 
2206  int n = (int) x;
2207 
2208  double tmp = x;
2209 
2210 
2211 
2212  while(--n){
2213 
2214  tmp -= 1.0;
2215 
2216  gamma *= tmp;
2217 
2218  }
2219 
2220  }
2221 
2222  else
2223 
2224  fprintf(stderr,"%lf is not of integer or half order\n",x);
2225 
2226  return gamma;
2227 
2228  }
2229 
2230 
2231 
2232  /**
2233 
2234  \brief Calculate the \a n zeros, \a z, of the Jacobi polynomial, i.e.
2235 
2236  \f$ P_n^{\alpha,\beta}(z) = 0 \f$
2237 
2238 
2239 
2240  This routine is only value for \f$( \alpha > -1, \beta > -1)\f$
2241 
2242  and uses polynomial deflation in a Newton iteration
2243 
2244  */
2245 
2246 
2247 
2248  static void Jacobz(const int n, double *z, const double alpha,
2249 
2250  const double beta){
2251 
2252  register int i,j,k;
2253 
2254  double dth = M_PI/(2.0*(double)n);
2255 
2256  double poly,pder,rlast=0.0;
2257 
2258  double sum,delr,r;
2259 
2260  double one = 1.0, two = 2.0;
2261 
2262 
2263 
2264  if(!n)
2265 
2266  return;
2267 
2268 
2269 
2270  for(k = 0; k < n; ++k){
2271 
2272  r = -cos((two*(double)k + one) * dth);
2273 
2274  if(k) r = 0.5*(r + rlast);
2275 
2276 
2277 
2278  for(j = 1; j < STOP; ++j){
2279 
2280  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
2281 
2282 
2283 
2284  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
2285 
2286 
2287 
2288  delr = -poly / (pder - sum * poly);
2289 
2290  r += delr;
2291 
2292  if( fabs(delr) < EPS ) break;
2293 
2294  }
2295 
2296  z[k] = r;
2297 
2298  rlast = r;
2299 
2300  }
2301 
2302  return;
2303 
2304  }
2305 
2306 
2307 
2308 
2309 
2310  /**
2311 
2312  \brief Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal
2313 
2314  matrix from the three term recurrence relationship.
2315 
2316 
2317 
2318  Set up a symmetric tridiagonal matrix
2319 
2320 
2321 
2322  \f$ \left [ \begin{array}{ccccc}
2323 
2324  a[0] & b[0] & & & \\
2325 
2326  b[0] & a[1] & b[1] & & \\
2327 
2328  0 & \ddots & \ddots & \ddots & \\
2329 
2330  & & \ddots & \ddots & b[n-2] \\
2331 
2332  & & & b[n-2] & a[n-1] \end{array} \right ] \f$
2333 
2334 
2335 
2336  Where the coefficients a[n], b[n] come from the recurrence relation
2337 
2338 
2339 
2340  \f$ b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) \f$
2341 
2342 
2343 
2344  where \f$ j=n+1\f$ and \f$p_j(z)\f$ are the Jacobi (normalized)
2345 
2346  orthogonal polynomials \f$ \alpha,\beta > -1\f$( integer values and
2347 
2348  halves). Since the polynomials are orthonormalized, the tridiagonal
2349 
2350  matrix is guaranteed to be symmetric. The eigenvalues of this
2351 
2352  matrix are the zeros of the Jacobi polynomial.
2353 
2354  */
2355 
2356 
2357 
2358  void JacZeros(const int n, double *a, double*b, const double alpha,
2359 
2360  const double beta){
2361 
2362 
2363 
2364  int i,j;
2365 
2366  RecCoeff(n,a,b,alpha,beta);
2367 
2368 
2369 
2370  double **z = new double*[n];
2371 
2372  for(i = 0; i < n; i++)
2373 
2374  {
2375 
2376  z[i] = new double[n];
2377 
2378  for(j = 0; j < n; j++)
2379 
2380  {
2381 
2382  z[i][j] = 0.0;
2383 
2384  }
2385 
2386  }
2387 
2388  for(i = 0; i < n; i++)
2389 
2390  {
2391 
2392  z[i][i] = 1.0;
2393 
2394  }
2395 
2396 
2397 
2398  // find eigenvalues and eigenvectors
2399 
2400  TriQL(n, a, b,z);
2401 
2402 
2403 
2404  return;
2405 
2406  }
2407 
2408 
2409 
2410  /**
2411 
2412  \brief The routine finds the recurrence coefficients \a a and
2413 
2414  \a b of the orthogonal polynomials
2415 
2416  */
2417 
2418  static void RecCoeff(const int n, double *a, double *b,const double alpha,
2419 
2420  const double beta){
2421 
2422 
2423 
2424  int i;
2425 
2426  double apb, apbi,a2b2;
2427 
2428 
2429 
2430  if(!n)
2431 
2432  return;
2433 
2434 
2435 
2436  // generate normalised terms
2437 
2438  apb = alpha + beta;
2439 
2440  apbi = 2.0 + apb;
2441 
2442 
2443 
2444  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
2445 
2446  a[0] = (beta-alpha)/apbi;
2447 
2448  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2449 
2450 
2451 
2452  a2b2 = beta*beta-alpha*alpha;
2453 
2454 
2455 
2456  for(i = 1; i < n-1; i++){
2457 
2458  apbi = 2.0*(i+1) + apb;
2459 
2460  a[i] = a2b2/((apbi-2.0)*apbi);
2461 
2462  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2463 
2464  ((apbi*apbi-1)*apbi*apbi));
2465 
2466  }
2467 
2468 
2469 
2470  apbi = 2.0*n + apb;
2471 
2472  a[n-1] = a2b2/((apbi-2.0)*apbi);
2473 
2474 
2475 
2476  }
2477 
2478 
2479 
2480 
2481 
2482  /** \brief QL algorithm for symmetric tridiagonal matrix
2483 
2484 
2485 
2486  This subroutine is a translation of an algol procedure,
2487 
2488  num. math. \b 12, 377-383(1968) by martin and wilkinson, as modified
2489 
2490  in num. math. \b 15, 450(1970) by dubrulle. Handbook for
2491 
2492  auto. comp., vol.ii-linear algebra, 241-248(1971). This is a
2493 
2494  modified version from numerical recipes.
2495 
2496 
2497 
2498  This subroutine finds the eigenvalues and first components of the
2499 
2500  eigenvectors of a symmetric tridiagonal matrix by the implicit QL
2501 
2502  method.
2503 
2504 
2505 
2506  on input:
2507 
2508  - n is the order of the matrix;
2509 
2510  - d contains the diagonal elements of the input matrix;
2511 
2512  - e contains the subdiagonal elements of the input matrix
2513 
2514  in its first n-2 positions.
2515 
2516  - z is the n by n identity matrix
2517 
2518 
2519 
2520  on output:
2521 
2522 
2523 
2524  - d contains the eigenvalues in ascending order.
2525 
2526  - e contains the weight values - modifications of the first component
2527 
2528  of normalised eigenvectors
2529 
2530  */
2531 
2532 
2533 
2534  static void TriQL(const int n, double *d,double *e, double **z){
2535 
2536  int m,l,iter,i,k;
2537 
2538  double s,r,p,g,f,dd,c,b;
2539 
2540 
2541 
2542  double MuZero = e[0];
2543 
2544 
2545 
2546  // Renumber the elements of e
2547 
2548  for(i = 0; i < n-1; i++)
2549 
2550  {
2551 
2552  e[i] = sqrt(e[i+1]);
2553 
2554  }
2555 
2556  e[n-1] = 0.0;
2557 
2558 
2559 
2560 
2561 
2562  for (l=0;l<n;l++) {
2563 
2564  iter=0;
2565 
2566  do {
2567 
2568  for (m=l;m<n-1;m++) {
2569 
2570  dd=fabs(d[m])+fabs(d[m+1]);
2571 
2572  if (fabs(e[m])+dd == dd) break;
2573 
2574  }
2575 
2576  if (m != l) {
2577 
2578  if (iter++ == STOP){
2579 
2580  fprintf(stderr,"triQL: Too many iterations in TQLI");
2581 
2582  exit(1);
2583 
2584  }
2585 
2586  g=(d[l+1]-d[l])/(2.0*e[l]);
2587 
2588  r=sqrt((g*g)+1.0);
2589 
2590  g=d[m]-d[l]+e[l]/(g+sign(r,g));
2591 
2592  s=c=1.0;
2593 
2594  p=0.0;
2595 
2596  for (i=m-1;i>=l;i--) {
2597 
2598  f=s*e[i];
2599 
2600  b=c*e[i];
2601 
2602  if (fabs(f) >= fabs(g)) {
2603 
2604  c=g/f;
2605 
2606  r=sqrt((c*c)+1.0);
2607 
2608  e[i+1]=f*r;
2609 
2610  c *= (s=1.0/r);
2611 
2612  } else {
2613 
2614  s=f/g;
2615 
2616  r=sqrt((s*s)+1.0);
2617 
2618  e[i+1]=g*r;
2619 
2620  s *= (c=1.0/r);
2621 
2622  }
2623 
2624  g=d[i+1]-p;
2625 
2626  r=(d[i]-g)*s+2.0*c*b;
2627 
2628  p=s*r;
2629 
2630  d[i+1]=g+p;
2631 
2632  g=c*r-b;
2633 
2634 
2635 
2636  // Calculate the eigenvectors
2637 
2638  for(k = 0; k < n; k++)
2639 
2640  {
2641 
2642  f = z[k][i+1];
2643 
2644  z[k][i+1] = s*z[k][i] + c*f;
2645 
2646  z[k][i] = c*z[k][i] - s*f;
2647 
2648  }
2649 
2650 
2651 
2652  }
2653 
2654  d[l]=d[l]-p;
2655 
2656  e[l]=g;
2657 
2658  e[m]=0.0;
2659 
2660  }
2661 
2662  } while (m != l);
2663 
2664  }
2665 
2666 
2667 
2668  // order eigenvalues
2669 
2670  // Since we only need the first component of the eigenvectors
2671 
2672  // to calcualte the weight, we only swap the first components
2673 
2674  for(i = 0; i < n-1; ++i){
2675 
2676  k = i;
2677 
2678  p = d[i];
2679 
2680  for(l = i+1; l < n; ++l)
2681 
2682  if (d[l] < p) {
2683 
2684  k = l;
2685 
2686  p = d[l];
2687 
2688  }
2689 
2690  d[k] = d[i];
2691 
2692  d[i] = p;
2693 
2694 
2695 
2696  double temp = z[0][k];
2697 
2698  z[0][k] = z[0][i];
2699 
2700  z[0][i] = temp;
2701 
2702  }
2703 
2704 
2705 
2706  // Calculate the weights
2707 
2708  for(i =0 ; i < n; i++)
2709 
2710  {
2711 
2712  e[i] = MuZero*z[0][i]*z[0][i];
2713 
2714  }
2715 
2716  }
2717 
2718 
2719 
2720  /**
2721 
2722  \brief Calcualtes the Jacobi-kronrod matrix by determining the
2723 
2724  \a a and \b coefficients.
2725 
2726 
2727 
2728  The first \a 3n+1 coefficients are already known
2729 
2730 
2731 
2732  For more information refer to:
2733 
2734  "Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"
2735 
2736  */
2737 
2738  void JKMatrix(int n, double *a, double *b)
2739 
2740  {
2741 
2742  int i,j,k,m;
2743 
2744  // Working storage
2745 
2746  int size = (int)floor(n/2.0)+2;
2747 
2748  double *s = new double[size];
2749 
2750  double *t = new double[size];
2751 
2752 
2753 
2754  // Initialize s and t to zero
2755 
2756  for(i = 0; i < size; i++)
2757 
2758  {
2759 
2760  s[i] = 0.0;
2761 
2762  t[i] = 0.0;
2763 
2764  }
2765 
2766 
2767 
2768  t[1] = b[n+1];
2769 
2770  for(m = 0; m <= n-2; m++)
2771 
2772  {
2773 
2774  double u = 0.0;
2775 
2776  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
2777 
2778  {
2779 
2780  int l = m-k;
2781 
2782  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
2783 
2784  s[k+1] = u;
2785 
2786  }
2787 
2788 
2789 
2790  // Swap the contents of s and t
2791 
2792  double *hold = s;
2793 
2794  s = t;
2795 
2796  t = hold;
2797 
2798  }
2799 
2800 
2801 
2802 
2803 
2804  for(j = (int)floor(n/2.0); j >= 0; j--)
2805 
2806  {
2807 
2808  s[j+1] = s[j];
2809 
2810  }
2811 
2812 
2813 
2814  for(m = n-1; m <= 2*n-3; m++)
2815 
2816  {
2817 
2818  double u = 0;
2819 
2820  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
2821 
2822  {
2823 
2824  int l = m-k;
2825 
2826  j = n-1-l;
2827 
2828  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
2829 
2830  s[j+1] = u;
2831 
2832  }
2833 
2834 
2835 
2836  if(m%2 == 0)
2837 
2838  {
2839 
2840  k = m/2;
2841 
2842  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
2843 
2844 
2845 
2846  }else
2847 
2848  {
2849 
2850  k = (m+1)/2;
2851 
2852  b[k+n+1] = s[j+1]/s[j+2];
2853 
2854  }
2855 
2856 
2857 
2858 
2859 
2860  // Swap the contents of s and t
2861 
2862  double *hold = s;
2863 
2864  s = t;
2865 
2866  t = hold;
2867 
2868  }
2869 
2870 
2871 
2872  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
2873 
2874 
2875 
2876  }
2877 
2878 
2879 
2880  /**
2881 
2882  \brief
2883 
2884  Given a weight function \f$w(t)\f$ through the first \a n+1
2885 
2886  coefficients \a a and \a b of its orthogonal polynomials
2887 
2888  this routine generates the first \a n recurrence coefficients for the orthogonal
2889 
2890  polynomials relative to the modified weight function \f$(t-z)w(t)\f$.
2891 
2892 
2893 
2894  The result will be placed in the array \a a0 and \a b0.
2895 
2896  */
2897 
2898 
2899 
2900  void chri1(int n, double* a, double* b, double* a0,
2901 
2902  double* b0,double z)
2903 
2904  {
2905 
2906 
2907 
2908  double q = ceil(3.0*n/2);
2909 
2910  int size = (int)q+1;
2911 
2912  if(size < n+1)
2913 
2914  {
2915 
2916  fprintf(stderr,"input arrays a and b are too short\n");
2917 
2918  }
2919 
2920  double* r = new double[n+1];
2921 
2922  r[0] = z - a[0];
2923 
2924  r[1] = z - a[1] - b[1]/r[0];
2925 
2926  a0[0] = a[1] + r[1] - r[0];
2927 
2928  b0[0] = -r[0]*b[0];
2929 
2930 
2931 
2932  if(n == 1)
2933 
2934  {
2935 
2936  return;
2937 
2938  }
2939 
2940  int k = 0;
2941 
2942  for(k = 1; k < n; k++)
2943 
2944  {
2945 
2946  r[k+1] = z - a[k+1] - b[k+1]/r[k];
2947 
2948  a0[k] = a[k+1] + r[k+1] - r[k];
2949 
2950  b0[k] = b[k] * r[k]/r[k-1];
2951 
2952  }
2953 
2954 
2955  }
2956 
2957  /**
2958 
2959  \brief
2960 
2961  Calcualte the bessel function of the first kind with complex double input y.
2962  Taken from Numerical Recipies in C
2963 
2964  Returns a complex double
2965  */
2966 
2967 
2968  std::complex<Nektar::NekDouble> ImagBesselComp(int n,std::complex<Nektar::NekDouble> y)
2969  {
2970  std::complex<Nektar::NekDouble> z (1.0,0.0);
2971  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
2972  std::complex<Nektar::NekDouble> zarg;
2973  Nektar::NekDouble tol = 1e-15;
2974  int maxit = 10000;
2975  int i = 1;
2976 
2977  zarg = -0.25*y*y;
2978 
2979  while (abs(z) > tol && i <= maxit){
2980  z = z*(1.0/i/(i+n)*zarg);
2981  if (abs(z) <= tol) break;
2982  zbes = zbes + z;
2983  i++;
2984  }
2985  zarg = 0.5*y;
2986  for (i=1;i<=n;i++){
2987  zbes = zbes*zarg;
2988  }
2989  return zbes;
2990 
2991  }
2992 
2993 
2994 
2995 } // end of namespace
2996 
2997 
2998 
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:972
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2418
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:274
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through.
Definition: Polylib.cpp:1386
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:1744
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:346
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
Definition: Polylib.cpp:1513
void zwlk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:746
void JacZeros(const int n, double *a, double *b, const double alpha, const double beta)
Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal.
Definition: Polylib.cpp:2358
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2900
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:2120
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:80
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:19
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:1806
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2738
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:1058
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:420
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:1684
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:1158
double NekDouble
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:198
NekDouble L
The namespace associated with the the Polylib library (Polylib introduction)
Definition: Polylib.cpp:31
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
Definition: Polylib.cpp:37
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2534
static void Jacobz(const int n, double *z, const double alpha, const double beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e.
Definition: Polylib.cpp:2248
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Radau-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:554
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
Definition: Polylib.cpp:1578
std::complex< Nektar::NekDouble > ImagBesselComp(int n, std::complex< Nektar::NekDouble > y)
Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Rec...
Definition: Polylib.cpp:2968
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
Definition: Polylib.cpp:1448
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:23
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:1260
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:1626
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:57
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:1920
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2174
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:142