Polylib.cpp 43.2 KB
 Blake Nelson committed Dec 17, 2009 1 2 3 4 5 6 7 8 9 10 11  #include #include #include #include "Polylib.h" #include  Mike committed Apr 24, 2015 12 13 14 15 16 #include #include  Blake Nelson committed Dec 17, 2009 17 18 19 20 21 22 /// Maximum number of iterations in polynomial defalation routine Jacobz #define STOP 30 /// Precision tolerance for two points to be similar  Hui Xu committed May 26, 2017 23 #define EPS 100*DBL_EPSILON  Blake Nelson committed Dec 17, 2009 24 25 26 27 28 29 30 31 32  /// return the sign(b)*a #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a)) namespace Polylib {  Hui Xu committed Apr 28, 2017 33 34 35 36  /// The following function is used to circumvent/reduce "Subtractive Cancellation" /// The expression 1/dz is replaced by optinvsub(.,.) /// Added on 26 April 2017  Hui Xu committed Jun 04, 2017 37 38  double optdiff(double xl, double xr) {  Hui Xu committed May 26, 2017 39  double m_xln, m_xrn;  Hui Xu committed Jun 04, 2017 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56  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); } }  Hui Xu committed Jun 04, 2017 57  double laginterp(double z, int j, const double *zj, int np)  Hui Xu committed Jun 04, 2017 58 59 60 61 62 63 64 65  { 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); 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) { 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]; }  Mike committed Apr 24, 2015 2954   Blake Nelson committed Dec 17, 2009 2955 2956  }  Mike committed Apr 24, 2015 2957 2958 2959 2960 2961 2962 2963 2964 2965 2966 2967 2968 2969 2970 2971 2972 2973 2974 2975 2976 2977 2978 2979 2980 2981 2982 2983 2984 2985 2986 2987 2988 2989 2990 2991 2992  /** \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; }  Blake Nelson committed Dec 17, 2009 2993 2994 2995 2996 2997 2998  } // end of namespace