Package com.polytechnik.utils
Class OrthogonalPolynomialsBasisFunctionsCalculatable<T extends BasisFunctionsCalculatable>
java.lang.Object
com.polytechnik.utils.OrthogonalPolynomialsBasisFunctionsCalculatable<T>
- Direct Known Subclasses:
OrthogonalPolynomialsChebyshevBasis
,OrthogonalPolynomialsHermiteEBasis
,OrthogonalPolynomialsLaguerreBasis
,OrthogonalPolynomialsLegendreBasis
,OrthogonalPolynomialsLegendreShiftedBasis
,OrthogonalPolynomialsMonomialsBasis
,OrthogonalPolynomialsRecurrenceABBasis
public class OrthogonalPolynomialsBasisFunctionsCalculatable<T extends BasisFunctionsCalculatable>
extends Object
Construct orthogonal polynomials in a given basis.
-
Field Summary
Fields -
Constructor Summary
Constructors -
Method Summary
Modifier and TypeMethodDescriptionprivate double
Dfx
(double[] m, double x) Calculate f'(x)=\sum_{k=0}^{k=n} m[k] Q'_k(x).private double
fx
(double[] m, double x) Calculate f(x)=\sum_{k=0}^{k=n} m[k] Q_k(x).getAB
(double[][] p) Calculate a and b recurrence coefficients from orthogonal polynomials.getABFromQuadrature
(double[] x, double[] w) Calculate a and b recurrence coefficients from a Gauss quadrature.getABRadauModified
(int n, RecurrenceAB ab_in) Modify b[n-2] to have Jacobi matrix generate (n-1)-th Radau polynomial exactly.double[]
getChristoffelOFunctionFromMatrix
(int n, double[] QQ) 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.double[]
getChristoffelOFunctionFromPolynomials
(int n, double[][] p) 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.protected EVSolver
Override it for a differentEVSolver
implementation (e.g.getGaussQuadratureFromAB
(int n, RecurrenceAB ab) Corresponding to quadratures built on n-th polynomial.getGaussQuadratureFromMoments
(int n, double[] moments) Calculate Gauss quadrature from moments in Qk(x) basis.getGaussQuadratureFromMoments_Direct
(int n, double[] moments) Similar togetGaussQuadratureFromMoments
but calculates quadratures directly: find n-th orthogonal polynomial, then find roots, and finally weights.getGaussQuadratureFromPolynomials
(int n, double[][] p) Corresponding to quadratures built on p[n] polynomial.static double[][]
getGramSchmidtOrthogonalized
(int n, double[] QQ) Given <Q_i|Q_j> matrix calculate first n orthonormal vectors using Gram Schmidt orthogonalization.getKronrodQuadrature
(int n, double[] moments) Calculates Gauss-Kronrod quadratures.getLebesgueQuadratureFromMatrices
(int n, double[] QQf, double[] QQ) Calculate Lebesgue quadrature.getLebesgueQuadratureFromMoments
(int n, double[] fmoments, double[] moments) Calculate Lebesgue quadrature from moments.double[]
Recover 0 ..double[]
getMomentsFromPolynomials
(double[][] p) Recover 0 ..double[]
getMultipleOrthogonalPolynomial_TypeII
(int[] multiindex, double[][] measures) Calculates single Type II multiple orthogonal polynomial.static double[][]
getOrthogonalized
(int n, double[] QQ) The same asgetGramSchmidtOrthogonalized
, slower, but hopefully? more stable.double[][]
getOrthogonalPolynomialsFirstKind
(int n, double[] moments) Calculate first n orhogonal polynomials (0..n-1) from first 2n-1 moments.double[][]
double[][]
getOrthogonalPolynomialsFromAB
(int kind, RecurrenceAB ab) Build orthogonal polynomials from a and b recurrence coefficients.double[][]
Polynomial root finder in original basis.getQuadraturesForMultipleOrthogonalPolynomial
(int[] multiindex, double[][] measures) Calculates Gaussian quadratures for multiple orthogonal polynimials.double[]
getRadauMomentsFromMoments
(double[] moments) Generate Radau <(x-x0)Q_k> moments from <Q_k> moments.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.double[][]
getRadauPolynomials
(double[][] p) getRadauQuadratureFromMoments
(int n, double[] moments) Corresponding to quadratures built on n+1 -th Radau polynomial.getRadauQuadratureFromPolynomials
(int n, double[][] p) Corresponding to quadratures built on Radau polynomial, constructed from p (measue: (x-x0)dmu).double
Calculate the root of getRadauMultiplied([1]).double[][]
getSecondKindOrthogonalPolynomialsFromFirstKind
(double[][] p) double
getSignToMatch
(int k) Orthogonal polynomials are defined within a factor of +-1 multiplied by.double[]
getSingleOrthogonalPolynomial
(int n, double[] moments) Calculates a single orthogonal polynomial by solving a linear system.double[]
getWeights
(double[] nodespolynomial, double[] x, double[][] p) Weights by integration.double[]
getWeightsByIntegration
(double[] nodespolynomial, double[] x, double[] moments) Weights by integration.double[]
getWeightsChristoffel
(double[] x, double[][] p) Weights using Christoffel function.getXNodes
(double[] nodespolynomial) Find the roots by solving polynomial equation nodespolynomial(X)=0.double[]
getXNodesFromAB
(int n, RecurrenceAB ab) Find roots by finding the eigenvalues of Jacobi matrix, n<ab.numAk().double[]
getxP
(double[] p) Multiply P(x)=\sum_{k=0}^{k=n} p[k] Q_k(x) by x in Q_k(x) basis.static void
printAB
(String header, RecurrenceAB ab) Used for diagnostics.static void
printQuadratures
(String txt, int n, double[] moments, OrthogonalPolynomialsBasisFunctionsCalculatable<? extends BasisFunctionsCalculatable> q) Used for diagnostics.
-
Field Details
-
B
Basis functions. -
EV
-
-
Constructor Details
-
OrthogonalPolynomialsBasisFunctionsCalculatable
- Parameters:
B
- The basis in which to construct orthogonal polynomials. We construct some other orthogonal polynomials in the basisB
.
-
-
Method Details
-
getEVSolver
Override it for a differentEVSolver
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 usinggetQQMatr
, inverts is, and callsgetKK
.- 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> bygetQQMatr()
usingsetQlQmExpansionCoefs()
.- 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 asgetGramSchmidtOrthogonalized
, slower, but hopefully? more stable. The method is equivalent to repeatedly callinggetSingleOrthogonalPolynomial
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
Calculate a and b recurrence coefficients from a Gauss quadrature. Used for unit test. -
getAB
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:
Usage example obtain three term recurrence coefficients from generalized moments (modified moments) inx*p_{k}=a_{k+1}*p_{k+1}+b_{k}p_{k}+a_{k}p_{k-1}
Chebyshev
basisvar 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
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
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
-
getOrthogonalPolynomialsFirstKindFromAB
-
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 usegetMomentsFromAB
. -
getMomentsFromAB
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
Find roots by finding the eigenvalues of Jacobi matrix, n<ab.numAk(). -
getXNodes
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 basisB
.- 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
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 togetGaussQuadratureFromMoments
but calculates quadratures directly: find n-th orthogonal polynomial, then find roots, and finally weights. Used for accuracy estimation ofgetGaussQuadratureFromMoments
.- 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
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 bygetAB()
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 bygetQQMatr()
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
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
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
Used for diagnostics.
-