Package com.polytechnik.utils
Class RecurrenceAB
java.lang.Object
com.polytechnik.utils.RecurrenceAB
- All Implemented Interfaces:
BasisFunctionsCalculatable
,BasisFunctionsMultipliable
,BasisPolynomials
,ConfederateMatrixCalculatable
,SimpleBasisPolynomials
,ThreeTermRecurrenceCalculatable
- Direct Known Subclasses:
RecurrenceABWithMultiplicationCached
public class RecurrenceAB
extends Object
implements BasisPolynomials, ThreeTermRecurrenceCalculatable
A normalized basis
constructed from given three--term recurrence coefficients.
The basis uses three--term recurrence coefficients obtained
analytically or from sampled data by
e.g
getAB()
method.
A typical usage: given three--term recurrence coefficients
obtain corresponding orthogonal polynomial
in other polynomials basis, see convertBasisToPBASIS
as an example
var ab=new RecurrenceAB(a,b);
double [][] polynomials_in_PBASIS=ab.convertBasisToPBASIS(ab.numAk(),PBASIS);
// now polynomials_in_PBASIS[k] is k-th orthogobal polynomial with (a,b) recurrence.
A common case of PBASIS
is
Monomials()
,
Chebyshev()
, etc.
Another example: assume we have two polynomial bases known
only through their three--term recurrence coefficients (a1,b1) and (a2,b2),
then we can transform
corresponding orthogonal polynomials through each other.
See such an example in the unit test TestOrthogonalPolynomialsBasisFunctionsCalculatable.test_AB_conversion_PBASIS()
method.
// create random recurrence coefficients a1 b1, a2 b2 from random measures.
var r=new java.util.Random(1234);
var n=10;
var MC=new OrthogonalPolynomialsLegendreBasis();
var tmp1=MC.getAB(MC.getOrthogonalPolynomialsFirstKind(n,MC.B.getMomentsFromSample(2*n-1,IntStream.range(0,1000).mapToDouble(o-> -1+2*r.nextDouble()).toArray())));
var a1=tmp1.getAllAk();
var b1=tmp1.getAllBk();
var tmp2=MC.getAB(MC.getOrthogonalPolynomialsFirstKind(n,MC.B.getMomentsFromSample(2*n-1,IntStream.range(0,1000).mapToDouble(o-> -1+2*r.nextDouble()).toArray())));
var a2=tmp2.getAllAk();
var b2=tmp2.getAllBk();
Then run the code:
RecurrenceAB ab1=new RecurrenceAB(a1,b1),ab2=new RecurrenceAB(a2,b2);
// the pols2_in_pols1Basis_direct[k] is the k-th orthogonal polynomial of ab2 in the basis of ab1.
double [][] pols2_in_pols1Basis_direct=ab2.convertBasisToPBASIS(ab2.numAk(),ab1);
// same polynomial can be obtained via an intermediate, e.g. the basis
var M=new OrthogonalPolynomialsChebyshevBasis();
double [][] pols1=M.getOrthogonalPolynomialsFirstKindFromAB(ab1), // ab1 polynomials in Chebyshev basis
pols2=M.getOrthogonalPolynomialsFirstKindFromAB(ab2); // ab2 polynomials in Chebyshev basis
for(int k=0;k<pols2.length;k++){
// this is the k-th orthogonal polynomial of ab2 in the basis of ab1.
double [] pols2k_in_pols1_via_intermediate=TMatr.expandOverLowerDiagBasis(pols2[k],pols1);
// it matches exactly the pols2_in_pols1Basis_direct[k].
System.err.println("direct="+java.util.Arrays.toString(pols2_in_pols1Basis_direct[k]));
System.err.println("via intermediate="+java.util.Arrays.toString(pols2k_in_pols1_via_intermediate));
}
-
Field Summary
Fields -
Constructor Summary
Constructors -
Method Summary
Modifier and TypeMethodDescriptiondouble
calculate
(int n, double x) Calculate \( Q_n(x) \) using three term recurrence at given x.double
getAk
(int k) a_kdouble[]
getAllAk()
Convenience method.double[]
getAllBk()
Convenience method.double[]
getaxbP
(double[] p, double ma, double mb) Do not mistake ma and mb multipliers with a and b recurrence coefficients.protected double[][]
getBasisFunctionsDifferentiated
(int numQl) Calculate the derivative of basis functions.double[][]
getBasisFunctionsMultipliedByPolynomial
(int numQl, double[] pol) Used as multiplication foundation when an analytic formula forsetQlQmExpansionCoefs(int, int, double[])
is not available, callsgetBasisFunctionsMultipliedByPolynomial(numQ,pol,this)
double[][]
getBasisFunctionsMultipliedByPolynomial
(int numQl, double[] pol, BasisFunctionsCalculatable PBASIS) Re-implement base interface method not to usePBASIS.mult2Pol
becasuse forRecurrenceAB
this method is the foundation of multiplication operation; instead, we now repeatedly apply therecurrence
formula to array.double[][]
getBasisFunctionsOnPolynomialArgument
(int numQl, double[] pol, BasisFunctionsCalculatable PBASIS) This method calculates \( Q_l(pol(x))\), where pol(x) is a given polynomial and \( Q_l \) is basis finction of this (notPBASIS
!) basis.double
getBk
(int k) b_kdouble[]
getConfederateMatrix
(double[] coefs) Confederate matrix used inPolynomialRootsConfederateMatrix
.double
getD0
(double[] psi, double x) Calculate psi_jQ_j(x) using three term recurrence and Clenshaw recurrence formula.double
getD1
(double[] psi, double x) First derivative at x.double
getD2
(double[] psi, double x) Second derivative at x.double[]
getDifferentiated
(double[] Tcoefs) Differentiate a function in basis,double[]
getIntegrated
(double[] Tcoefs) Calculate the integral as inverse to differentiation operation.double[]
getMomentsFromSample
(int nmoments, double[] x, double[] w, int nelsinsample) Given observations x and weights w calculate the 0..nmoms-1 moments in the basis.double[]
getNext
(int n, double[] Tnm1, double[] Tnm2, boolean flag_type2, BasisFunctionsCalculatable PBASIS) A mockup non-implemented method fromBasisPolynomials
, it is put there with other non-optimized default methods fromBasisPolynomials
.double[]
getNext
(int n, double[] Tnm1, double[] Tnm2, BasisFunctionsCalculatable PBASIS) Obtain next polynomial in PBASIS according to three term recurrence coefficients.double[]
getPQkMomentsFromQkMoments
(double[] P, double[] moments) Calculate <P(x) Q_k(x)> moments from <Q_k(x)> moments.static void
Unit test a general 3 term recurrence basis.double[]
mult2Pol
(double[] p1, double[] p2) Multiply two polynomials in basis usingBasisFunctionsMultipliable.setQlQmExpansionCoefs(int,int,double [])
operator.int
numAk()
The number of a_k available.sdiv1
(double[] p, double x0) Synthetically divide p by (x-x0).void
setQlQmExpansionCoefs
(int l, int m, double[] s) CallsgetBasisFunctionsMultipliedByPolynomial
, theBasisPolynomials#setQlQmExpansionCoefs
typically uses analytic expression; this implementation instead uses three term recurrence implemented ingetBasisFunctionsMultipliedByPolynomial
.Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
Methods inherited from interface com.polytechnik.utils.BasisFunctionsCalculatable
convertFunctionTo_Basis_from_Monomials, convertFunctionTo_Basis_from_PBASIS, convertFunctionTo_Monomials_from_Basis, getMomentsFromSample, getMomentsFromSample, getONE, getX0Moments, getXQkMomentsFromQkMoments
Methods inherited from interface com.polytechnik.utils.BasisFunctionsMultipliable
getKK, getMomentsOfMeasureProducingPolynomialInKK, getMomentsOfMeasureProducingPolynomialInKK_MQQM, getQQMatr, getTwoQuadraticFormsProductAsQuadraticForm, sdiv
Methods inherited from interface com.polytechnik.utils.BasisPolynomials
setNewtonBinomialLikeCoefs
Methods inherited from interface com.polytechnik.utils.SimpleBasisPolynomials
convertBasisToPBASIS, convertFunctionTo_PBASIS_from_Basis, getBasisFunctionsOnPolynomialArgument
-
Field Details
-
a
private final double[] a -
b
private final double[] b
-
-
Constructor Details
-
RecurrenceAB
public RecurrenceAB(double[] a, double[] b)
-
-
Method Details
-
getAk
-
getBk
-
numAk
public int numAk()The number of a_k available.- Specified by:
numAk
in interfaceThreeTermRecurrenceCalculatable
-
getAllAk
public double[] getAllAk()Convenience method. -
getAllBk
public double[] getAllBk()Convenience method. -
calculate
public double calculate(int n, double x) Calculate \( Q_n(x) \) using three term recurrence at given x.- Specified by:
calculate
in interfaceBasisFunctionsCalculatable
- Parameters:
n
- Polynomial index (starts with 0).x
- The point to calculate the polynomial.- Returns:
- The value of the polynomial Tn at x.
-
getD0
public double getD0(double[] psi, double x) Calculate psi_jQ_j(x) using three term recurrence and Clenshaw recurrence formula. See L. Fox, I. B. Parker (1968), Chebyshev Polynomials in Numerical Analysis, page 56.- Specified by:
getD0
in interfaceBasisFunctionsCalculatable
- Parameters:
psi
- Coefficients.x
- Point to estimate the sum.- Returns:
- The sum.
-
getD1
public double getD1(double[] psi, double x) First derivative at x.- Specified by:
getD1
in interfaceBasisFunctionsCalculatable
- Parameters:
psi
- Coefficients.x
- Point to estimate the sum.- Returns:
- The derivative.
-
getD2
public double getD2(double[] psi, double x) Second derivative at x.- Specified by:
getD2
in interfaceBasisFunctionsCalculatable
- Parameters:
psi
- Coefficients.x
- Point to estimate the sum.- Returns:
- The second derivative.
-
getBasisFunctionsDifferentiated
protected double[][] getBasisFunctionsDifferentiated(int numQl) Calculate the derivative of basis functions. Used bygetDifferentiated
andgetIntegrated
.- Parameters:
numQl
- The number of basis functions to calculate the derivative.- Returns:
- An array of Q'_0, Q'_1,... Q'_{numQl-1}. For an implementation with internal caching the method may return more functions than the numQl.
-
getDifferentiated
public double[] getDifferentiated(double[] Tcoefs) Description copied from interface:BasisFunctionsCalculatable
Differentiate a function in basis,- Specified by:
getDifferentiated
in interfaceBasisFunctionsCalculatable
- Parameters:
Tcoefs
- : \( F(x)=\sum_k Q_k(x) Qcoefs[k] \).- Returns:
- \( dF/dx=\sum_k Q'_k(x) Qcoefs[k] \).
-
getIntegrated
public double[] getIntegrated(double[] Tcoefs) Calculate the integral as inverse to differentiation operation.- Specified by:
getIntegrated
in interfaceBasisFunctionsCalculatable
- Parameters:
Tcoefs
- : \( F(x)=\sum_k Q_k(x) Qcoefs[k] \).- Returns:
- \( \int F(x)dx=\sum_k Qcoefs[k] \int Q_k(x) dx \). The integral is typically normalized as I(0)=0.
-
getaxbP
public double[] getaxbP(double[] p, double ma, double mb) Do not mistake ma and mb multipliers with a and b recurrence coefficients.- Specified by:
getaxbP
in interfaceBasisFunctionsCalculatable
- Parameters:
p
- Input polynomial.- Returns:
- p*(ma*x+mb).
-
getNext
Obtain next polynomial in PBASIS according to three term recurrence coefficients.- Specified by:
getNext
in interfaceBasisPolynomials
- Specified by:
getNext
in interfaceSimpleBasisPolynomials
- Parameters:
n
- Q index.Tnm1
- Input \( Q_{n-1} \).Tnm2
- Input \( Q_{n-2} \).PBASIS
- The basis in which to calculte.- Returns:
- \( Q_{n} \) polynomial of this basis in PBASIS basis, array dimension n+1.
-
getNext
public double[] getNext(int n, double[] Tnm1, double[] Tnm2, boolean flag_type2, BasisFunctionsCalculatable PBASIS) A mockup non-implemented method fromBasisPolynomials
, it is put there with other non-optimized default methods fromBasisPolynomials
. For type 2 polynomials it is recommended to create a new instance ofRecurrenceAB(double [] a,double [] b)
witha
andb
arrays having0
-th element removed. If the arguments used to create an instance are not availble -- one can get the values from convenience methodsgetAllAk()
andgetAllBk()
.- Specified by:
getNext
in interfaceBasisPolynomials
- Parameters:
n
- T index.Tnm1
- Input \( Q_{n-1} \).Tnm2
- Input \( Q_{n-2} \).flag_type2
- Whether to use type 2 polynomials. If type 2 polynomials are not available --- throw an exception. If not implemented -- an exception will be thrown. This method is typically not the proper way to obtain type 2 polynomials, use thegetOrthogonalPolynomialsSecondKindFromAB
instead.PBASIS
- The basis in which to calculte, form this basis we only need P*(a*x+b) operation, that is performed relatively this bais recurrence.- Returns:
- \( Q_{n} \) polynomial of this basis in PBASIS basis, array dimension n+1.
-
getBasisFunctionsOnPolynomialArgument
public double[][] getBasisFunctionsOnPolynomialArgument(int numQl, double[] pol, BasisFunctionsCalculatable PBASIS) Description copied from interface:SimpleBasisPolynomials
This method calculates \( Q_l(pol(x))\), where pol(x) is a given polynomial and \( Q_l \) is basis finction of this (notPBASIS
!) basis. This method may not be stable for a large dimensions ofpol.length
. Equivalent to numQl applications ofSimpleBasisPolynomials.getNext(int, double[], double[], com.polytechnik.utils.BasisFunctionsCalculatable)
topol
array. The method can be used to calculate e.g. \( T_n(P_k(x)) \) or \( T_n(P_k(L_m(x^j))) \). For example
meansnew Legendre().getBasisFunctionsOnPolynomialArgument(4,new double []{0,0,0,0,1},new Chebyshev()) ={ { 1.0 }, { 0.0, 0.0, 0.0, 0.0, 1.0 }, { 0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.75 }, { 0.0, 0.0, 0.0, 0.0, 0.375, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.625 } }
\( P_0(T_4(x))=T_0 \)
\( P_1(T_4(x))=T_4(x) \)
\( P_2(T_4(x))=0.25*T_0+0.75*T_8(x) \)
\( P_3(T_4(x))=0.375*T_4(x)+0.625*T_{12}(x) \).- Specified by:
getBasisFunctionsOnPolynomialArgument
in interfaceSimpleBasisPolynomials
- Parameters:
numQl
- The number of basis functions to calculate.pol
- The polynomial to use in \( Q_l(pol(x))\), it is given in thePBASIS
basis.PBASIS
- the basis in whichpol
is given and the result is returned.- Returns:
- An array of numQl dimension, elements: \( Q_0(pol(x)) \),
\( Q_1(pol(x)) \), ... \( Q_{numQl-1}(pol(x)) \) in the basis
PBASIS
.
-
getConfederateMatrix
public double[] getConfederateMatrix(double[] coefs) Confederate matrix used inPolynomialRootsConfederateMatrix
. Corresponds to \[ p_{k}=(x/a_{k}-b_{k-1}/a_{k})*p_{k-1} - (a_{k-1}/a_{k})*p_{k-2} \] recurrence.- Specified by:
getConfederateMatrix
in interfaceConfederateMatrixCalculatable
- Parameters:
coefs
- Polynomial coefficients.- Returns:
- confederate matrix.
-
setQlQmExpansionCoefs
public void setQlQmExpansionCoefs(int l, int m, double[] s) CallsgetBasisFunctionsMultipliedByPolynomial
, theBasisPolynomials#setQlQmExpansionCoefs
typically uses analytic expression; this implementation instead uses three term recurrence implemented ingetBasisFunctionsMultipliedByPolynomial
.- Specified by:
setQlQmExpansionCoefs
in interfaceBasisFunctionsMultipliable
- Parameters:
l
- l.m
- m.s
- A buffer of dimension at leastl+m+1
to save the expansion coeficients.
-
getBasisFunctionsMultipliedByPolynomial
public double[][] getBasisFunctionsMultipliedByPolynomial(int numQl, double[] pol, BasisFunctionsCalculatable PBASIS) Re-implement base interface method not to usePBASIS.mult2Pol
becasuse forRecurrenceAB
this method is the foundation of multiplication operation; instead, we now repeatedly apply therecurrence
formula to array.- Specified by:
getBasisFunctionsMultipliedByPolynomial
in interfaceSimpleBasisPolynomials
- Parameters:
numQl
- The number ofQl*pol
to calculate.pol
- Polynomial to multiply, it is given in basisPBASIS
.PBASIS
- The basis in whichpol
is given and the result is returned.- Returns:
- An array of
numQl
dimension, elements:Q0*pol, Q1*pol, ... Q_{numQl-1}*pol
of this basis functionsQl
multiplied bypol
; the result is returned in basisPBASIS
.
-
getBasisFunctionsMultipliedByPolynomial
public double[][] getBasisFunctionsMultipliedByPolynomial(int numQl, double[] pol) Used as multiplication foundation when an analytic formula forsetQlQmExpansionCoefs(int, int, double[])
is not available, callsgetBasisFunctionsMultipliedByPolynomial(numQ,pol,this)
- Specified by:
getBasisFunctionsMultipliedByPolynomial
in interfaceBasisFunctionsMultipliable
- Parameters:
numQl
- The number of Ql*P to calculate.pol
- Polynomial to multiply.- Returns:
- An array of numQl dimension, elements: Q0*pol, Q1*pol, ... Q_{numQl-1}*pol. For example
new Chebyshev().getBasisFunctionsMultipliedByPolynomial(3,new double[]{0,0,0,0,0,1}) ={ { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 }, { 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5 }, { 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5 } }
-
getPQkMomentsFromQkMoments
public double[] getPQkMomentsFromQkMoments(double[] P, double[] moments) Description copied from interface:BasisFunctionsMultipliable
Calculate <P(x) Q_k(x)> moments from <Q_k(x)> moments.- Specified by:
getPQkMomentsFromQkMoments
in interfaceBasisFunctionsMultipliable
- Parameters:
moments
- <F|Qk>.- Returns:
- <F|P(x)Qk>.
-
mult2Pol
public double[] mult2Pol(double[] p1, double[] p2) Description copied from interface:BasisFunctionsMultipliable
Multiply two polynomials in basis usingBasisFunctionsMultipliable.setQlQmExpansionCoefs(int,int,double [])
operator.- Specified by:
mult2Pol
in interfaceBasisFunctionsMultipliable
- Parameters:
p1
- Firts polynomial.p2
- Second polynomial.- Returns:
- p*q. For example:
.new Chebyshev().mult2Pol(new double[]{0,0,1},new double[]{0,0,0,0,0,1}) ={ 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5 }
-
sdiv1
Description copied from interface:BasisFunctionsCalculatable
Synthetically divide p by (x-x0). Calculations in basis.- Specified by:
sdiv1
in interfaceBasisFunctionsCalculatable
- Parameters:
p
- Input polynomial;x0
- To divide on (x-x0)- Returns:
- Quotient q and remainder r in class SDivRes.
-
getMomentsFromSample
public double[] getMomentsFromSample(int nmoments, double[] x, double[] w, int nelsinsample) Description copied from interface:BasisFunctionsCalculatable
Given observations x and weights w calculate the 0..nmoms-1 moments in the basis.- Specified by:
getMomentsFromSample
in interfaceBasisFunctionsCalculatable
- Parameters:
nmoments
- The number of moments to calculate.x
- observation samples.w
- weights.nelsinsample
- How many elements to sample.- Returns:
- nelsinsample generalzied moments, \( k=0\dots nmoms-1 \): \( \sum_{j=0}^{nelsinsample-1} Q_k(x[j])*w[j] \).
-
main
Unit test a general 3 term recurrence basis.
-