#include #include #include #include "Polylib.h" #include #include #include /// Maximum number of iterations in polynomial defalation routine Jacobz #define STOP 30 /// Precision tolerance for two points to be similar #define EPS 100*DBL_EPSILON /// return the sign(b)*a #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a)) namespace Polylib { /// The following function is used to circumvent/reduce "Subtractive Cancellation" /// The expression 1/dz is replaced by optinvsub(.,.) /// Added on 26 April 2017 double optdiff(double xl, double xr) { double m_xln, m_xrn; int m_expn; int m_digits = static_cast(fabs(floor(log10(DBL_EPSILON)))-1); if (fabs(xl-xr)<1.e-4){ m_expn = static_cast(floor(log10(fabs(xl-xr)))); m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn)); return powl(10.0L,-m_digits)*(m_xln-m_xrn); }else{ return (xl-xr); } } double laginterp(double z, int j, const double *zj, int np) { double temp = 1.0; for (int i=0; i -1, \beta > -1 \f$and its derivative at the \a np points in \a z[i] - If \a poly_in = NULL then only calculate derivatice - If \a polyd = NULL then only calculate polynomial - To calculate the polynomial this routine uses the recursion relationship (see appendix A ref [4]) : \f$ \begin{array}{rcl} P^{\alpha,\beta}_0(z) &=& 1 \\ P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\ a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z) P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\ a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\ a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\ a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1) (2n + \alpha + \beta + 2) \\ a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2) \end{array} \f$- To calculate the derivative of the polynomial this routine uses the relationship (see appendix A ref [4]) : \f$ \begin{array}{rcl} b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z) + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\ b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\ b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\ b^3_n(z) &=& 2(n+\alpha)(n+\beta) \end{array} \f$- Note the derivative from this routine is only valid for -1 < \a z < 1. */ void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta){ register int i; double zero = 0.0, one = 1.0, two = 2.0; if(!np) return; if(n == 0){ if(poly_in) for(i = 0; i < np; ++i) poly_in[i] = one; if(polyd) for(i = 0; i < np; ++i) polyd[i] = zero; } else if (n == 1){ if(poly_in) for(i = 0; i < np; ++i) poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]); if(polyd) for(i = 0; i < np; ++i) polyd[i] = 0.5*(alpha + beta + two); } else{ register int k; double a1,a2,a3,a4; double two = 2.0, apb = alpha + beta; double *poly, *polyn1,*polyn2; if(poly_in){ // switch for case of no poynomial function return polyn1 = (double *)malloc(2*np*sizeof(double)); polyn2 = polyn1+np; poly = poly_in; } else{ polyn1 = (double *)malloc(3*np*sizeof(double)); polyn2 = polyn1+np; poly = polyn2+np; } for(i = 0; i < np; ++i){ polyn2[i] = one; polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]); } for(k = 2; k <= n; ++k){ a1 = two*k*(k + apb)*(two*k + apb - two); a2 = (two*k + apb - one)*(alpha*alpha - beta*beta); a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb); a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb); a2 /= a1; a3 /= a1; a4 /= a1; for(i = 0; i < np; ++i){ poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i]; polyn2[i] = polyn1[i]; polyn1[i] = poly [i]; } } if(polyd){ a1 = n*(alpha - beta); a2 = n*(two*n + alpha + beta); a3 = two*(n + alpha)*(n + beta); a4 = (two*n + alpha + beta); a1 /= a4; a2 /= a4; a3 /= a4; // note polyn2 points to polyn1 at end of poly iterations for(i = 0; i < np; ++i){ polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i]; polyd[i] /= (one - z[i]*z[i]); } } free(polyn1); } return; } /** \brief Calculate the derivative of Jacobi polynomials \li Generates a vector \a poly of values of the derivative of the \a n th order Jacobi polynomial \f$ P^(\alpha,\beta)_n(z)\f$at the \a np points \a z. \li To do this we have used the relation \n \f$ \frac{d}{dz} P^{\alpha,\beta}_n(z) = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) \f$\li This formulation is valid for \f$ -1 \leq z \leq 1 \f$*/ void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta) { register int i; double one = 1.0; if(n == 0) for(i = 0; i < np; ++i) polyd[i] = 0.0; else{ //jacobf(np,z,polyd,n-1,alpha+one,beta+one); jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one); for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one); } return; } /** \brief Calculate the Gamma function , \f$ \Gamma(n)\f$, for integer values and halves. Determine the value of \f$\Gamma(n)\f$using: \f$ \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$where \f$ \Gamma(1/2) = \sqrt(\pi)\f$*/ double gammaF(const double x){ double gamma = 1.0; if (x == -0.5) gamma = -2.0*sqrt(M_PI); else if (!x) return gamma; else if ((x-(int)x) == 0.5){ int n = (int) x; double tmp = x; gamma = sqrt(M_PI); while(n--){ tmp -= 1.0; gamma *= tmp; } } else if ((x-(int)x) == 0.0){ int n = (int) x; double tmp = x; while(--n){ tmp -= 1.0; gamma *= tmp; } } else fprintf(stderr,"%lf is not of integer or half order\n",x); return gamma; } /** \brief Calculate the \a n zeros, \a z, of the Jacobi polynomial, i.e. \f$ P_n^{\alpha,\beta}(z) = 0 \f$This routine is only value for \f$( \alpha > -1, \beta > -1)\f$and uses polynomial deflation in a Newton iteration */ static void Jacobz(const int n, double *z, const double alpha, const double beta){ register int i,j,k; double dth = M_PI/(2.0*(double)n); double poly,pder,rlast=0.0; double sum,delr,r; double one = 1.0, two = 2.0; if(!n) return; for(k = 0; k < n; ++k){ r = -cos((two*(double)k + one) * dth); if(k) r = 0.5*(r + rlast); for(j = 1; j < STOP; ++j){ jacobfd(1,&r,&poly, &pder, n, alpha, beta); for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]); delr = -poly / (pder - sum * poly); r += delr; if( fabs(delr) < EPS ) break; } z[k] = r; rlast = r; } return; } /** \brief Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from the three term recurrence relationship. Set up a symmetric tridiagonal matrix \f$ \left [ \begin{array}{ccccc} a[0] & b[0] & & & \\ b[0] & a[1] & b[1] & & \\ 0 & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & b[n-2] \\ & & & b[n-2] & a[n-1] \end{array} \right ] \f$Where the coefficients a[n], b[n] come from the recurrence relation \f$ b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) \f$where \f$ j=n+1\f$and \f$p_j(z)\f$are the Jacobi (normalized) orthogonal polynomials \f$ \alpha,\beta > -1\f$( integer values and halves). Since the polynomials are orthonormalized, the tridiagonal matrix is guaranteed to be symmetric. The eigenvalues of this matrix are the zeros of the Jacobi polynomial. */ void JacZeros(const int n, double *a, double*b, const double alpha, const double beta){ int i,j; RecCoeff(n,a,b,alpha,beta); double **z = new double*[n]; for(i = 0; i < n; i++) { z[i] = new double[n]; for(j = 0; j < n; j++) { z[i][j] = 0.0; } } for(i = 0; i < n; i++) { z[i][i] = 1.0; } // find eigenvalues and eigenvectors TriQL(n, a, b,z); delete[] z; return; } /** \brief The routine finds the recurrence coefficients \a a and \a b of the orthogonal polynomials */ static void RecCoeff(const int n, double *a, double *b,const double alpha, const double beta){ int i; double apb, apbi,a2b2; if(!n) return; // generate normalised terms apb = alpha + beta; apbi = 2.0 + apb; b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero a[0] = (beta-alpha)/apbi; b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi)); a2b2 = beta*beta-alpha*alpha; for(i = 1; i < n-1; i++){ apbi = 2.0*(i+1) + apb; a[i] = a2b2/((apbi-2.0)*apbi); b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/ ((apbi*apbi-1)*apbi*apbi)); } apbi = 2.0*n + apb; a[n-1] = a2b2/((apbi-2.0)*apbi); } /** \brief QL algorithm for symmetric tridiagonal matrix This subroutine is a translation of an algol procedure, num. math. \b 12, 377-383(1968) by martin and wilkinson, as modified in num. math. \b 15, 450(1970) by dubrulle. Handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). This is a modified version from numerical recipes. This subroutine finds the eigenvalues and first components of the eigenvectors of a symmetric tridiagonal matrix by the implicit QL method. on input: - n is the order of the matrix; - d contains the diagonal elements of the input matrix; - e contains the subdiagonal elements of the input matrix in its first n-2 positions. - z is the n by n identity matrix on output: - d contains the eigenvalues in ascending order. - e contains the weight values - modifications of the first component of normalised eigenvectors */ static void TriQL(const int n, double *d,double *e, double **z){ int m,l,iter,i,k; double s,r,p,g,f,dd,c,b; double MuZero = e[0]; // Renumber the elements of e for(i = 0; i < n-1; i++) { e[i] = sqrt(e[i+1]); } e[n-1] = 0.0; for (l=0;l=l;i--) { f=s*e[i]; b=c*e[i]; if (fabs(f) >= fabs(g)) { c=g/f; r=sqrt((c*c)+1.0); e[i+1]=f*r; c *= (s=1.0/r); } else { s=f/g; r=sqrt((s*s)+1.0); e[i+1]=g*r; s *= (c=1.0/r); } g=d[i+1]-p; r=(d[i]-g)*s+2.0*c*b; p=s*r; d[i+1]=g+p; g=c*r-b; // Calculate the eigenvectors for(k = 0; k < n; k++) { f = z[k][i+1]; z[k][i+1] = s*z[k][i] + c*f; z[k][i] = c*z[k][i] - s*f; } } d[l]=d[l]-p; e[l]=g; e[m]=0.0; } } while (m != l); } // order eigenvalues // Since we only need the first component of the eigenvectors // to calcualte the weight, we only swap the first components for(i = 0; i < n-1; ++i){ k = i; p = d[i]; for(l = i+1; l < n; ++l) if (d[l] < p) { k = l; p = d[l]; } d[k] = d[i]; d[i] = p; double temp = z[0][k]; z[0][k] = z[0][i]; z[0][i] = temp; } // Calculate the weights for(i =0 ; i < n; i++) { e[i] = MuZero*z[0][i]*z[0][i]; } } /** \brief Calcualtes the Jacobi-kronrod matrix by determining the \a a and \b coefficients. The first \a 3n+1 coefficients are already known For more information refer to: "Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules" */ void JKMatrix(int n, double *a, double *b) { int i,j,k,m; // Working storage int size = (int)floor(n/2.0)+2; double *s = new double[size]; double *t = new double[size]; // Initialize s and t to zero for(i = 0; i < size; i++) { s[i] = 0.0; t[i] = 0.0; } t[1] = b[n+1]; for(m = 0; m <= n-2; m++) { double u = 0.0; for(k = (int)floor((m+1)/2.0); k >= 0; k--) { int l = m-k; u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1]; s[k+1] = u; } // Swap the contents of s and t double *hold = s; s = t; t = hold; } for(j = (int)floor(n/2.0); j >= 0; j--) { s[j+1] = s[j]; } for(m = n-1; m <= 2*n-3; m++) { double u = 0; for(k = m+1-n; k <= floor((m-1)/2.0); k++) { int l = m-k; j = n-1-l; u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2]; s[j+1] = u; } if(m%2 == 0) { k = m/2; a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2]; }else { k = (m+1)/2; b[k+n+1] = s[j+1]/s[j+2]; } // Swap the contents of s and t double *hold = s; s = t; t = hold; } a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1]; } /** \brief Given a weight function \f$w(t)\f$through the first \a n+1 coefficients \a a and \a b of its orthogonal polynomials this routine generates the first \a n recurrence coefficients for the orthogonal polynomials relative to the modified weight function \f$(t-z)w(t)\f\$. The result will be placed in the array \a a0 and \a b0. */ void chri1(int n, double* a, double* b, double* a0, double* b0,double z) { double q = ceil(3.0*n/2); int size = (int)q+1; if(size < n+1) { fprintf(stderr,"input arrays a and b are too short\n"); } double* r = new double[n+1]; r[0] = z - a[0]; r[1] = z - a[1] - b[1]/r[0]; a0[0] = a[1] + r[1] - r[0]; b0[0] = -r[0]*b[0]; if(n == 1) { delete[] r; return; } int k = 0; for(k = 1; k < n; k++) { r[k+1] = z - a[k+1] - b[k+1]/r[k]; a0[k] = a[k+1] + r[k+1] - r[k]; b0[k] = b[k] * r[k]/r[k-1]; } delete[] r; } /** \brief Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C Returns a complex double */ std::complex ImagBesselComp(int n,std::complex y) { std::complex z (1.0,0.0); std::complex zbes (1.0,0.0); std::complex zarg; Nektar::NekDouble tol = 1e-15; int maxit = 10000; int i = 1; zarg = -0.25*y*y; while (abs(z) > tol && i <= maxit){ z = z*(1.0/i/(i+n)*zarg); if (abs(z) <= tol) break; zbes = zbes + z; i++; } zarg = 0.5*y; for (i=1;i<=n;i++){ zbes = zbes*zarg; } return zbes; } } // end of namespace