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 Details

    • a

      private final double[] a
    • b

      private final double[] b
  • Constructor Details

    • RecurrenceAB

      public RecurrenceAB(double[] a, double[] b)
  • Method Details

    • getAk

      public double getAk(int k)
      a_k
      Specified by:
      getAk in interface ThreeTermRecurrenceCalculatable
    • getBk

      public double getBk(int k)
      b_k
      Specified by:
      getBk in interface ThreeTermRecurrenceCalculatable
    • numAk

      public int numAk()
      The number of a_k available.
      Specified by:
      numAk in interface ThreeTermRecurrenceCalculatable
    • 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 interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      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 by getDifferentiated and getIntegrated.
      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 interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      Parameters:
      p - Input polynomial.
      Returns:
      p*(ma*x+mb).
    • getNext

      public double[] getNext(int n, double[] Tnm1, double[] Tnm2, BasisFunctionsCalculatable PBASIS)
      Obtain next polynomial in PBASIS according to three term recurrence coefficients.
      Specified by:
      getNext in interface BasisPolynomials
      Specified by:
      getNext in interface SimpleBasisPolynomials
      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 from BasisPolynomials, it is put there with other non-optimized default methods from BasisPolynomials. For type 2 polynomials it is recommended to create a new instance of RecurrenceAB(double [] a,double [] b) with a and b arrays having 0-th element removed. If the arguments used to create an instance are not availble -- one can get the values from convenience methods getAllAk() and getAllBk().
      Specified by:
      getNext in interface BasisPolynomials
      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 the getOrthogonalPolynomialsSecondKindFromAB 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 (not PBASIS!) basis. This method may not be stable for a large dimensions of pol.length. Equivalent to numQl applications of SimpleBasisPolynomials.getNext(int, double[], double[], com.polytechnik.utils.BasisFunctionsCalculatable) to pol 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
      
      new 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 }
      }
      
      means
      \( 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 interface SimpleBasisPolynomials
      Parameters:
      numQl - The number of basis functions to calculate.
      pol - The polynomial to use in \( Q_l(pol(x))\), it is given in the PBASIS basis.
      PBASIS - the basis in which pol 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 in PolynomialRootsConfederateMatrix. 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 interface ConfederateMatrixCalculatable
      Parameters:
      coefs - Polynomial coefficients.
      Returns:
      confederate matrix.
    • setQlQmExpansionCoefs

      public void setQlQmExpansionCoefs(int l, int m, double[] s)
      Calls getBasisFunctionsMultipliedByPolynomial, the BasisPolynomials#setQlQmExpansionCoefs typically uses analytic expression; this implementation instead uses three term recurrence implemented in getBasisFunctionsMultipliedByPolynomial.
      Specified by:
      setQlQmExpansionCoefs in interface BasisFunctionsMultipliable
      Parameters:
      l - l.
      m - m.
      s - A buffer of dimension at least l+m+1 to save the expansion coeficients.
    • getBasisFunctionsMultipliedByPolynomial

      public double[][] getBasisFunctionsMultipliedByPolynomial(int numQl, double[] pol, BasisFunctionsCalculatable PBASIS)
      Re-implement base interface method not to use PBASIS.mult2Pol becasuse for RecurrenceAB this method is the foundation of multiplication operation; instead, we now repeatedly apply the recurrence formula to array.
      Specified by:
      getBasisFunctionsMultipliedByPolynomial in interface SimpleBasisPolynomials
      Parameters:
      numQl - The number of Ql*pol to calculate.
      pol - Polynomial to multiply, it is given in basis PBASIS.
      PBASIS - The basis in which pol 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 functions Ql multiplied by pol; the result is returned in basis PBASIS.
    • getBasisFunctionsMultipliedByPolynomial

      public double[][] getBasisFunctionsMultipliedByPolynomial(int numQl, double[] pol)
      Used as multiplication foundation when an analytic formula for setQlQmExpansionCoefs(int, int, double[]) is not available, calls getBasisFunctionsMultipliedByPolynomial(numQ,pol,this)
      Specified by:
      getBasisFunctionsMultipliedByPolynomial in interface BasisFunctionsMultipliable
      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 interface BasisFunctionsMultipliable
      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 using BasisFunctionsMultipliable.setQlQmExpansionCoefs(int,int,double []) operator.
      Specified by:
      mult2Pol in interface BasisFunctionsMultipliable
      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

      public SDivRes sdiv1(double[] p, double x0)
      Description copied from interface: BasisFunctionsCalculatable
      Synthetically divide p by (x-x0). Calculations in basis.
      Specified by:
      sdiv1 in interface BasisFunctionsCalculatable
      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 interface BasisFunctionsCalculatable
      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

      public static void main(String[] args)
      Unit test a general 3 term recurrence basis.