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