From: Martin Kronbichler Date: Fri, 6 Apr 2018 09:06:06 +0000 (+0200) Subject: Add user-visible implementation of Jacobi polynomial. X-Git-Tag: v9.0.0-rc1~195^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d4cbeaa7a669af764ec0e2f8732c1cc04fd68426;p=dealii.git Add user-visible implementation of Jacobi polynomial. --- diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 129eb33050..9246a20b48 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -709,6 +709,43 @@ namespace Polynomials static std::vector > generate_complete_basis (const unsigned int degree); }; + + + + /* + * Evaluate a Jacobi polynomial $ P_n^{\alpha, \beta}(x) $ specified by the + * parameters @p alpha, @p beta, @p n, where @p n is the degree of the + * Jacobi polynomial. + * + * @note The Jacobi polynomials are not orthonormal and are defined on the + * unit interval $[0, 1]$ as usual for deal.II, rather than $[-1, +1]$ often + * used in literature. @p x is the point of evaluation. + */ + template + Number + jacobi_polynomial_value(const unsigned int degree, + const int alpha, + const int beta, + const Number x); + + + /** + * Compute the roots of the Jacobi polynomials on the unit interval $[0, 1]$ + * of the given degree. These roots are used in several places inside the + * deal.II library, such as the Gauss-Lobatto quadrature formula or for the + * Hermite-like interpolation. + * + * The algorithm uses a Newton algorithm, using the zeros of the Chebyshev + * polynomials as an initial guess. This code has been tested for alpha and + * beta equal to zero (Legendre case), one (Gauss-Lobatto case) as well as + * two, so be careful when using it for other values as the Newton iteration + * might or might not converge. + */ + template + std::vector + jacobi_polynomial_roots(const unsigned int degree, + const int alpha, + const int beta); } @@ -790,6 +827,126 @@ namespace Polynomials ar &lagrange_weight; } + + + template + Number + jacobi_polynomial_value(const unsigned int degree, + const int alpha, + const int beta, + const Number x) + { + Assert(alpha >=0 && beta >= 0, + ExcNotImplemented("Negative alpha/beta coefficients not supported")); + // the Jacobi polynomial is evaluated using a recursion formula. + Number p0, p1; + + // The recursion formula is defined for the interval [-1, 1], so rescale + // to that interval here + const Number xeval = Number(-1) + 2. * x; + + // initial values P_0(x), P_1(x): + p0 = 1.0; + if (degree==0) + return p0; + p1 = ((alpha+beta+2)*xeval + (alpha-beta))/2; + if (degree==1) + return p1; + + for (unsigned int i=1; i + std::vector + jacobi_polynomial_roots(const unsigned int degree, + const int alpha, + const int beta) + { + std::vector x(degree, 0.5); + + // compute zeros with a Newton algorithm. + + // Set tolerance. For long double we might not always get the additional + // precision in a run time environment (e.g. with valgrind), so we must + // limit the tolerance to double. Since we do a Newton iteration, doing + // one more iteration after the residual has indicated convergence will be + // enough for all number types due to the quadratic convergence of + // Newton's method + + const Number tolerance + = 4 * std::max(static_cast(std::numeric_limits::epsilon()), + std::numeric_limits::epsilon()); + + // The following implementation follows closely the one given in the + // appendix of the book by Karniadakis and Sherwin: Spectral/hp element + // methods for computational fluid dynamics (Oxford University Press, + // 2005) + + // If symmetric, we only need to compute the half of points + const unsigned int n_points = (alpha == beta ? degree/2 : degree); + for (unsigned int k=0; k (2*k+1)/(2*degree) * + numbers::PI ); + if (k>0) + r = (r + x[k-1])/2; + + unsigned int converged = numbers::invalid_unsigned_int; + for (unsigned int i=1; i<1000; ++i) + { + Number s = 0.; + for (unsigned int i=0; i