Class OrthogonalPolynomialsBasisFunctionsCalculatable<T extends BasisFunctionsCalculatable>

java.lang.Object
com.polytechnik.utils.OrthogonalPolynomialsBasisFunctionsCalculatable<T>
  • Field Details

  • Constructor Details

    • OrthogonalPolynomialsBasisFunctionsCalculatable

      public OrthogonalPolynomialsBasisFunctionsCalculatable(T B)
      Parameters:
      B - The basis in which to construct orthogonal polynomials. We construct some other orthogonal polynomials in the basis B.
  • Method Details

    • getEVSolver

      protected EVSolver getEVSolver()
      Override it for a different EVSolver implementation (e.g. EVSolver_JNI_lapacke, direct JNI call to http://www.netlib.org/lapack/lapacke.html ) because my com/polytechnik/lapack may have problems.
    • getDensityMatrixProducingGivenPolynomial

      public LebesgueQuadratureWithEVData getDensityMatrixProducingGivenPolynomial(int n, double[] P_toproduce, double[] QQ)
      For a given polynomial P_toproduce generate an operator (density matrix), such that P_toproduce(x)=sum_i lambda^{[i]} (psi^{[i]}(x))^2 Such an approach allows to use Spur(f|rho) instead of <f|P_toproduce>, matrix spur allows further generalization, @see https://arxiv.org/abs/1807.06007.
      Parameters:
      n - the dimension.
      P_toproduce - a polynomial to produce, must have 0..2n-2, total 2n-1 array elements.
      QQ - scalar product matrix.
      Returns:
      Lebsgue quadrature, eigenvalues/eigenvectors of which produce given polynomial.
    • getChristoffelOFunctionFromMoments

      public double[] getChristoffelOFunctionFromMoments(int n, double[] moments)
      Calculates Christoffel function 1/K(x) as K(x)=Q(x)(<QQ>)^{-1}Q(x) polynomial from <Qk> moments. The method constructs QQ using getQQMatr, inverts is, and calls getKK.
      Parameters:
      n - the dimension.
      moments - The <Qk> moments, dim>=2*n-1.
      Returns:
      The polynomial Q(x)(<QQ>)^{-1}Q(x).
    • getChristoffelOFunctionFromMatrix

      public double[] getChristoffelOFunctionFromMatrix(int n, double[] QQ)
    • getChristoffelOFunctionFromPolynomials

      public double[] getChristoffelOFunctionFromPolynomials(int n, double[][] p)
    • getOrthogonalPolynomialsFirstKind

      public double[][] getOrthogonalPolynomialsFirstKind(int n, double[] moments)
      Calculate first n orhogonal polynomials (0..n-1) from first 2n-1 moments. All we need is <Q_l Q_m> matrix which is internally constructed from mom[k]=<Q_k> by getQQMatr() using setQlQmExpansionCoefs().
      Parameters:
      n - The number of polynomials to generate.
      moments - Moments, must be at least 2n-1 elements (0..2n-2 moments).
      Returns:
      Orthogonal Polynomials.
    • getSingleOrthogonalPolynomial

      public double[] getSingleOrthogonalPolynomial(int n, double[] moments)
      Calculates a single orthogonal polynomial by solving a linear system. Useful for unit tests.
      Parameters:
      n - Order of the polynomial.
      moments - Moments in basis, should have at least 2*n elements in array (0..2n-1 moments).
      Returns:
      Polynomial, an array of n+1 length, it is normalized as P[n]=1.
    • getOrthogonalized

      public static double[][] getOrthogonalized(int n, double[] QQ)
      The same as getGramSchmidtOrthogonalized, slower, but hopefully? more stable. The method is equivalent to repeatedly calling getSingleOrthogonalPolynomial with subsequent normalizing. Calculations can be performed in various basises Q_k.
      Parameters:
      n - The number of basis vectors to build.
      QQ - The <Q_i Q_j> matrix, i,j=0..n-1, dim=n*n.
    • getSignToMatch

      public double getSignToMatch(int k)
      Orthogonal polynomials are defined within a factor of +-1 multiplied by. This method return this factor to match standard definition, corresponding positive sign of the coefficient with maximal x^n, when alfa_n Q_n(x) is converted to monomials.
      Returns:
      +1 or -1.
    • getABFromQuadrature

      public RecurrenceAB getABFromQuadrature(double[] x, double[] w)
      Calculate a and b recurrence coefficients from a Gauss quadrature. Used for unit test.
    • getAB

      public RecurrenceAB getAB(double[][] p)
      Calculate a and b recurrence coefficients from orthogonal polynomials. The definitions as in Walter Gautschi, "Orthogonal Polynomials:Computation and approximation", page 12. Here a[k]=sqrt(beta_k) and b[k]=alfa_k. And the recurrence is:
      x*p_{k}=a_{k+1}*p_{k+1}+b_{k}p_{k}+a_{k}p_{k-1}
      Usage example obtain three term recurrence coefficients from generalized moments (modified moments) in Chebyshev basis
      
      var M=new OrthogonalPolynomialsChebyshevBasis();
      var ab=M.getAB(M.getOrthogonalPolynomialsFirstKind(7,new double []{1,0,0,0,0,0,0,0,0,0,0,0,0}))
      ab.getAllAk() ={ 1.0, 0.7071067811865476, 0.5, 0.5, 0.5, 0.5, 0.5 }
      ab.getAllBk() ={ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
      // first kind
      M.getOrthogonalPolynomialsFirstKindFromAB(ab)
      ={ { 1.0 },
         { 0.0, 1.414213562373095 },
         { 0.0, 0.0, 1.414213562373095 },
         { 0.0, 0.0, 0.0, 1.414213562373095 },
         { 0.0, 0.0, 0.0, 0.0, 1.414213562373095 },
         { 0.0, 0.0, 0.0, 0.0, 0.0, 1.414213562373095 },
         { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.414213562373095 }
      }
      // second kind
      var p2=M.getOrthogonalPolynomialsSecondKindFromAB(ab)
      ={ { 1.414213562373095 },
         { 0.0, 2.82842712474619 },
         { 1.414213562373095, 0.0, 2.82842712474619 },
         { 0.0, 2.82842712474619, 0.0, 2.82842712474619 },
         { 1.414213562373095, 0.0, 2.82842712474619, 0.0, 2.82842712474619 },
         { 0.0, 2.82842712474619, 0.0, 2.82842712474619, 0.0, 2.82842712474619 }
      }
      var ab2=M.getAB(p2)
      ab2.getAllAk()={ 0.7071067811865476, 0.5, 0.5, 0.5, 0.5, 0.5 }
      ab2.getAllBk()={ 0.0, 0.0, 0.0, 0.0, 0.0 }
       
      Parameters:
      p - Orthogonal polynomials.
      Returns:
      Three--term recurrece coefficients.
    • getABRadauModified

      public RecurrenceAB getABRadauModified(int n, RecurrenceAB ab_in)
      Modify b[n-2] to have Jacobi matrix generate (n-1)-th Radau polynomial exactly. See "Journal of Computational and applied Mathematics", Volume 127, Numbers 1|2, 15 January 2001, D.P. Laurie, "Computation of Gauss-type quadrature formulas", p.201-217.
    • getOrthogonalPolynomialsFromAB

      public double[][] getOrthogonalPolynomialsFromAB(int kind, RecurrenceAB ab)
      Build orthogonal polynomials from a and b recurrence coefficients. May not always be stable, but useful for second kind polynomials generation.
      Parameters:
      kind - , 0-first kind, 1-second kind, etc.
      ab - recurrence coefficients.
      Returns:
      An array of orthogonal polynomials, the number of polynomials calculated is ab.numAk()-kind.
    • getOrthogonalPolynomialsSecondKindFromAB

      public double[][] getOrthogonalPolynomialsSecondKindFromAB(RecurrenceAB ab)
    • getOrthogonalPolynomialsFirstKindFromAB

      public double[][] getOrthogonalPolynomialsFirstKindFromAB(RecurrenceAB ab)
    • getSecondKindOrthogonalPolynomialsFromFirstKind

      public double[][] getSecondKindOrthogonalPolynomialsFromFirstKind(double[][] p)
    • getMomentsFromPolynomials

      public double[] getMomentsFromPolynomials(double[][] p)
      Recover 0 .. 2n-1 moments from 0 .. n orthogonal polynomials. The 2n-th moment is the normalizing factor of p[n] and cannot be recovered. Used for unit tests. For production use getMomentsFromAB.
    • getMomentsFromAB

      public double[] getMomentsFromAB(RecurrenceAB ab)
      Recover 0 .. 2n-1 moments from a[0..n] b[0..n-1]. Used for unit tests and to recover moments for second kind orthogonal polynomials measure.
    • getWeights

      public double[] getWeights(double[] nodespolynomial, double[] x, double[][] p)
      Weights by integration. Integration is replaced by synthetic division by basis[] polynomials. The coefficient by p[0] give the value of integral. The weights are normalized as sum(weights)=<1>. Require all orthogonal polynomials to be available.
    • getWeightsChristoffel

      public double[] getWeightsChristoffel(double[] x, double[][] p)
      Weights using Christoffel function. Here x are the roots of p[x.length], and the weights are normalized as sum(weights)=<1>. Require two orthogonal polynomials to be available: p[n] & p[n-1], and p[0] for normalizing.
    • getWeightsByIntegration

      public double[] getWeightsByIntegration(double[] nodespolynomial, double[] x, double[] moments)
      Weights by integration. Here the x[] are the roots of nodespolynomial polynomial. Synthetically divide nodespolynomial by (x-x[k]), then integrate Legendre interpolatory polynomial directly using known values of moments[k]=<Q_k(x)>. Typically stable (except for exponentially small weights). The weights are normalized as sum(weights)=<1>.
      Parameters:
      nodespolynomial - Nodes polynomial.
      x - Points.
      moments - Moments, at least nodespolynomial.length-1 elements must present.
      Returns:
      The weights.
    • getXNodesFromAB

      public double[] getXNodesFromAB(int n, RecurrenceAB ab)
      Find roots by finding the eigenvalues of Jacobi matrix, n<ab.numAk().
    • getXNodes

      public PolynomialRootsFinder.PRoots getXNodes(double[] nodespolynomial)
      Find the roots by solving polynomial equation nodespolynomial(X)=0. The method calls the roots finded, then sorts the nodes according to real part of the root. The result of PolynomialRootsFinder may be unsorted.
      Parameters:
      nodespolynomial - A polynomial in basis B.
      Returns:
      The result of getXNodes: all roots are found (otherwise NaN) and then are sorted.
    • getGaussQuadratureFromPolynomials

      public Quadrature.GaussQuadratureWithNodesPolynomialData getGaussQuadratureFromPolynomials(int n, double[][] p)
      Corresponding to quadratures built on p[n] polynomial. The weights are normalized as sum(weights)=<1>. Used for unit tests.
    • getRadauPolynomials

      public double[][] getRadauPolynomials(double[][] p)
    • getRadauQuadratureFromPolynomials

      public Quadrature.GaussQuadratureWithNodesPolynomialData getRadauQuadratureFromPolynomials(int n, double[][] p)
      Corresponding to quadratures built on Radau polynomial, constructed from p (measue: (x-x0)dmu). The weights are normalized as sum(weights)=<1>. This approach may be less stable than direct calculation from moments. but it allows us to calculate quadratures for second kind polynomials.
      Parameters:
      n - The order of generated Radau quadrature is n+1.
      p - Orthogonal polynomials, there are should be at least n+2 polynomials.
      Returns:
      nodes weighs and nodes polynomial
    • getRadauMomentsFromMoments

      public double[] getRadauMomentsFromMoments(double[] moments)
      Generate Radau <(x-x0)Q_k> moments from <Q_k> moments.
    • getLebesgueQuadratureFromMatrices

      public LebesgueQuadratureWithEVData getLebesgueQuadratureFromMatrices(int n, double[] QQf, double[] QQ)
      Calculate Lebesgue quadrature. If QQf=QQX the result is Gaussian quadrature, otherwise it can be considered as Lebesgue quadratures.
      Parameters:
      n - Points number.
      QQf - <fQjQk> matrix.
      QQ - <QjQk> matrix.
      Returns:
      value-nodes, weights, eigenvaluies/eigenvectors.
    • getLebesgueQuadratureFromMoments

      public LebesgueQuadratureWithEVData getLebesgueQuadratureFromMoments(int n, double[] fmoments, double[] moments)
      Calculate Lebesgue quadrature from moments. If fmoments=Xmoments the result is Gaussian quadrature, otherwise it can be considered as Lebesgue quadratures.
      Parameters:
      n - Points number.
      fmoments - <fQj> moments.
      moments - <Qj> moemnts.
      Returns:
      value-nodes (eigenvalues), weights, eigenvaluies/eigenvectors array,QQf,QQ.
    • getGaussQuadratureFromMoments

      public Quadrature.GaussQuadratureWithEVData getGaussQuadratureFromMoments(int n, double[] moments)
      Calculate Gauss quadrature from moments in Qk(x) basis. Weights are normalized as sum(weights)=<1>. Uses generalized eigenvalues expansion. Most stable.
      Parameters:
      n - Order of the quadrature.
      moments - Moments in basis, the array must have at least 2n elements (0..2n-1 moments).
      Returns:
      nodes, weights, eigenvaluies/eigenvectors array,XQQ,QQ.
    • getGaussQuadratureFromMoments_Direct

      public Quadrature.GaussQuadratureWithNodesPolynomialData getGaussQuadratureFromMoments_Direct(int n, double[] moments)
      Similar to getGaussQuadratureFromMoments but calculates quadratures directly: find n-th orthogonal polynomial, then find roots, and finally weights. Used for accuracy estimation of getGaussQuadratureFromMoments.
      Parameters:
      n - Ponts in the quadrature.
      moments - Moments in basis, must be at least 2n elements (0..2n-1 moments).
      Returns:
      nodes, weights, and nodes polynomial.
    • getGaussQuadratureFromAB

      public Quadrature getGaussQuadratureFromAB(int n, RecurrenceAB ab)
      Corresponding to quadratures built on n-th polynomial. Weights are normalized as sum(weights)=<1>. Uses eigenvectors expansion in orthogonal polynomials basis. Most stable.
      Parameters:
      n - Quadrature order.
      ab - Recurrence coefficient, must have at least n+1 element, e.g. obtained by getAB() method.
      Returns:
      The quadrature.
    • getRadauQuadratureFromMoments

      public Quadrature.GaussQuadratureWithNodesPolynomialData getRadauQuadratureFromMoments(int n, double[] moments)
      Corresponding to quadratures built on n+1 -th Radau polynomial. Weights are normalized as sum(weights)=<1>. Uses generalized eigenvalues expansion. Most stable.
      Parameters:
      n - The order of generated Radau quadrature is n+1.
      moments - Moments, should be at least 2n+1 elements.
      Returns:
      x-nodes, weights, nodes polynomial.
    • getGramSchmidtOrthogonalized

      public static double[][] getGramSchmidtOrthogonalized(int n, double[] QQ)
      Given <Q_i|Q_j> matrix calculate first n orthonormal vectors using Gram Schmidt orthogonalization. In polynomial basis gives first n orthogonal polynomials. Note: may be numerically unstable.
      Parameters:
      n - The number of basis vectors to generate.
      QQ - The <Q_i|Q_j> matrix, generated by getQQMatr() method from moments, i,j=0..n-1, dim=n*n.
      Returns:
      Othogonal basis obtained by applying Gram Schmidt orthogonalization to basis vectors in original order.
    • getMultipleOrthogonalPolynomial_TypeII

      public double[] getMultipleOrthogonalPolynomial_TypeII(int[] multiindex, double[][] measures)
      Calculates single Type II multiple orthogonal polynomial.
      Parameters:
      multiindex - Multiindex. Typically all elements are the same.
      measures - Moment measures <Q_k>_mu_s. Each s-th measure should have at least sum(multiindex)+multiindex[s] elements.
      Returns:
      Type II multiple orthogonal polynomial.
    • getQuadraturesForMultipleOrthogonalPolynomial

      public Quadrature.MultipleOrthogonalityQuadrature getQuadraturesForMultipleOrthogonalPolynomial(int[] multiindex, double[][] measures)
      Calculates Gaussian quadratures for multiple orthogonal polynimials.
      Parameters:
      multiindex - Multiindex, typically all elements are the same.
      measures - Moment measures <Q_k>_mu_s. Each s-th measure should have at least sum(multiindex)+multiindex[s] elements.
      Returns:
      Quadratures. of x nodes. imag(x), all equal to 0 on success, complex roots indicate failure to find proper quadrature. weight[s], s=0..multiindex.length-1; multiindex.length-1+3-th element: the nodes polynomial in basis. The results tested to match exactly the Carlos F. Borges,Numer. Math. 67: 271-288 (1994), On a class of Gauss-like quadrature rules. However the results are not stable for arbitrary measures, the method TestOrthogonalPolynomialsBasisFunctionsCalculatable.test_MultipleOrthogonality() fails for some randomly generated measures.
    • getKronrodQuadrature

      public Quadrature.KronrodQuadrature getKronrodQuadrature(int n, double[] moments)
      Calculates Gauss-Kronrod quadratures.
      Parameters:
      n - Order or original Gauss Quadrature, final quadrature has 2n+1 order.
      moments - Moments, must have at least 3n+2 elements (0..3n+1 moments).
      Returns:
      Quadratures. The 6 elements of Quadrature.KronrodQuadrature returned: Gauss x nodes. imag(x) of Gauss nodes, all equal to 0 on success. weights for Gauss nodes. Kronrod x nodes. imag(x) of Kronrod nodes, all equal to 0 on success. weights for Kronrod nodes.
    • getxP

      public double[] getxP(double[] p)
      Multiply P(x)=\sum_{k=0}^{k=n} p[k] Q_k(x) by x in Q_k(x) basis.
    • fx

      private double fx(double[] m, double x)
      Calculate f(x)=\sum_{k=0}^{k=n} m[k] Q_k(x).
    • Dfx

      private double Dfx(double[] m, double x)
      Calculate f'(x)=\sum_{k=0}^{k=n} m[k] Q'_k(x).
    • getPolynomialRootsFinderInBasis

      public PolynomialRootsFinder getPolynomialRootsFinderInBasis()
      Polynomial root finder in original basis. Used for unit testing the quadratures.
    • getRadauMultiplied

      public double[] getRadauMultiplied(double[] nodes0)
      Multiply nodes0 polynomial, P(x)=\sum_{k=0}^{k=n} nodes0[k] Q_k(x), by (x-x0) or (x0-x) , whatever is positive to obtain Radau nodes polynomial. The method has to be implemented if Radau quadrature is required, the default version throws an exception.
    • getRadauRoot

      public double getRadauRoot()
      Calculate the root of getRadauMultiplied([1]).
    • printQuadratures

      public static void printQuadratures(String txt, int n, double[] moments, OrthogonalPolynomialsBasisFunctionsCalculatable<? extends BasisFunctionsCalculatable> q)
      Used for diagnostics.
    • printAB

      public static void printAB(String header, RecurrenceAB ab)
      Used for diagnostics.