001package kernels; 002 003import java.util.Arrays; 004 005import org.apache.commons.math3.linear.BlockRealMatrix; 006import org.apache.commons.math3.linear.EigenDecomposition; 007import org.apache.commons.math3.linear.RealMatrix; 008import org.apache.commons.math3.linear.RealVector; 009 010/** 011 * An integration kernel \( K(\theta,\eta) \) corresponding to an operator on the one-particle space 012 * of an integrable theory. 013 * 014 * <p>This class is abstract and doesn't provide the actual kernel function; 015 * this is left to subclasses. Rather, this class provides conversion routines 016 * from the kernel function to matrix elements of the operator, and 017 * for approximately finding its spectral data.</p> 018 * 019 * <p>To that end, the following orthonormal system of piecewise constant bumps is used as a basis: 020 * \[ \varphi_j (x) = \begin{cases} 021 * h^{-1/2} \quad & \text{if } -R + j h < x < -R+(j+1)h , \\ 022 * 0 & \text{otherwise} , 023 * \end{cases} \] 024 * 025 * with \(j=0,\ldots,N-1\), where \(R>0\) is a given "cutoff", \(N\) the number of basis elements, and \(h=2R/N\). 026 * All approximation results for the eigenvectors, the spectrum, etc. of the operator will depend on these 027 * approximation parameters \(R\) and \(N\). </p> 028 * 029 */ 030public abstract class OneParticleKernel 031{ 032 033 private boolean evCheckEnabled = true; 034 035 /** 036 * Evaluates the kernel \( K(\theta,\eta) \) at a given position \( \theta,\eta \). 037 * 038 * @param theta first argument of the kernel 039 * @param eta second argument of the kernel 040 * @return value of the kernel at that position. 041 */ 042 public abstract double kernelValue(double theta, double eta); 043 044 /** 045 * Disable the plausibility check of decaying eigenfunctions for this kernel. 046 * This should be called once, when for this specific kernel it is not expected that the eigenfunctions decay. 047 * The check performed by {@link #checkPlausibility} will then be disabled for all method calls. 048 * 049 * @see #findLowestEigenvector 050 * @see #findSpectrum 051 */ 052 public void disableEVCheck() { 053 this.evCheckEnabled = false; 054 } 055 056 /** 057 * Retrieve the matrix element of this kernel between two basis vectors. 058 * 059 * The matrix elements \( (\varphi_i , A \varphi_j) \) are approximated using the midpoint rule. 060 * 061 * @param i index of the first basis vector 062 * @param j index of the second basis vector 063 * @param n number of basis vectors 064 * @param cutoff the "cutoff" \(R\) 065 * @return the matrix element \( (\varphi_i , A \varphi_j) \) 066 */ 067 protected double matrixElement(int i, int j, int n, double cutoff) { 068 double step = 2.0*cutoff/n; 069 double theta = -cutoff + step*(i + 0.5); 070 double eta = -cutoff + step*(j + 0.5); 071 double mtxe = step * kernelValue(theta, eta); 072 return mtxe; 073 } 074 075 076 /** 077 * Compute a matrix approximation of this operator, in the basis described above. 078 * 079 * @param n number of basis vectors 080 * @param cutoff the "cutoff" \(R\) 081 * @return a matrix approximation of this operator in the given basis 082 */ 083 public RealMatrix asRealMatrix(int n, double cutoff) { 084 RealMatrix matrix = new BlockRealMatrix(n, n); 085 for (int i = 0; i < n; i++) { 086 for (int j = 0; j <= i; j++) { 087 double mtxe = matrixElement(i, j, n, cutoff); 088 matrix.setEntry(i, j, mtxe); 089 matrix.setEntry(j, i, mtxe); 090 } 091 } 092 return matrix; 093 } 094 095 /** 096 * Compute an approximate eigen-decomposition of this operator, in the basis described above. 097 * 098 * @param n number of basis vectors 099 * @param cutoff the "cutoff" \(R\) 100 * @return eigen-decomposition of this operator in the given basis 101 */ 102 public EigenDecomposition findEigendec(int n, double cutoff) { 103 RealMatrix a = this.asRealMatrix(n, cutoff); 104 return new EigenDecomposition(a); 105 } 106 107 /** 108 * Check whether an eigen-decomposition of this operator is plausible, 109 * namely, whether its lowest eigenfunction decays at the end of the interval \([-R,R]\). 110 * 111 * This is useful for our main use case, when we want to find the eigenfunction for the lowest eigenvector 112 * independent of the cutoff \(R\). 113 * The method also checks whether the eigenvalues are all real. 114 * (We expect to deal only with symmetric operators.) 115 * 116 * @see #disableEVCheck() 117 * @see #findEigendec(int, double) 118 * 119 * @param ed the eigen-decomposition to verify 120 * @param n number of basis vectors 121 * @throws RuntimeException if the matrix has complex eigenvalues 122 * @throws RuntimeException if the lowest eigenfunction does not decay 123 */ 124 protected void checkPlausibility(EigenDecomposition ed, int n) { 125 if (ed.hasComplexEigenvalues()) { 126 throw new RuntimeException("complex eigenvalues detected"); 127 } 128 if (evCheckEnabled) { 129 RealVector v = ed.getEigenvector(n-1); 130 double eps = 0.0; 131 for (int i = 0; i < 5; i++) { 132 eps += Math.abs(v.getEntry(i)); 133 eps += Math.abs(v.getEntry(n-1-i)); 134 } 135 if (eps > 1e-4) { 136 throw new RuntimeException("lowest eigenfunction does not decay"); 137 } 138 } 139 } 140 141 /** 142 * Approximate the eigenfunction for the lowest eigenvalue of this operator, in the basis described above. 143 * 144 * @param n number of basis vectors 145 * @param cutoff the "cutoff" \(R\) 146 * @return the lowest eigenvector in the given basis 147 */ 148 public RealVector findLowestEigenvector(int n, double cutoff) { 149 RealMatrix a = this.asRealMatrix(n, cutoff); 150 EigenDecomposition e = new EigenDecomposition(a); 151 checkPlausibility(e, n); 152 return e.getEigenvector(a.getRowDimension()-1); 153 } 154 155 /** 156 * Approximate the spectrum of this operator, in the basis described above. 157 * The eigenvalues are returned in ascending order. 158 * 159 * @param n number of basis vectors 160 * @param cutoff the "cutoff" \(R\) 161 * @return the eigenvalues of the operator in this approximation 162 */ 163 public double[] findSpectrum(int n, double cutoff) { 164 RealMatrix a = this.asRealMatrix(n, cutoff); 165 EigenDecomposition e = new EigenDecomposition(a); 166 checkPlausibility(e, n); 167 double[] evs = e.getRealEigenvalues(); 168 Arrays.sort(evs); 169 return evs; 170 } 171 172 /** 173 * Approximate the lowest eigenvalue of this operator, using the basis described above. 174 * 175 * @param n number of basis vectors 176 * @param cutoff the "cutoff" \(R\) 177 * @return the lowest eigenvalue of the operator in this approximation 178 */ 179 public double findLowestEigenvalue(int n, double cutoff) { 180 return findSpectrum(n, cutoff)[0]; 181 } 182 183}