Nektar++
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
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
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
Definition: Polylib.cpp:932
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
Definition: Polylib.cpp:2444
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
Definition: Polylib.cpp:234
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:1340
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Definition: Polylib.cpp:1770
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
Definition: Polylib.cpp:306
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:1500
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:706
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:2384
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
Definition: Polylib.cpp:2926
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:2146
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:40
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:14
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points.
Definition: Polylib.cpp:1832
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
Definition: Polylib.cpp:2764
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
Definition: Polylib.cpp:1018
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:380
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Definition: Polylib.cpp:1710
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
Definition: Polylib.cpp:1116
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
Definition: Polylib.cpp:158
The namespace associated with the the Polylib library (Polylib introduction)
Definition: Polylib.cpp:26
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:2560
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:2274
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Definition: Polylib.cpp:514
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:1584
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:1416
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:18
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
Definition: Polylib.cpp:1216
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an.
Definition: Polylib.cpp:1652
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1946
double gammaF(const double)
Calculate the Gamma function , , for integer.
Definition: Polylib.cpp:2200
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:102