4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 1// Eigenvalues of a symmetric tridiagonal matrix (diagonal d, off-diagonal e,
2// e[i] couples i and i+1) by Sturm-sequence bisection. Zero off-diagonals
3// split it into independent blocks, which the count handles.
5function countBelow(d: Float64Array, e: Float64Array, x: number): number {
6 let count = 0;
7 let q = 1;
8 for (let i = 0; i < d.length; i++) {
9 const e2 = i > 0 ? e[i - 1] * e[i - 1] : 0;
10 q = d[i] - x - (i > 0 ? e2 / q : 0);
11 if (q === 0) q = -1e-300;
12 if (q < 0) count++;
13 }
14 return count;
15}
17export function tridiagEigenvalues(d: Float64Array, e: Float64Array, count = d.length): number[] {
18 const n = d.length;
19 if (n === 0) return [];
20 let lo = Infinity;
21 let hi = -Infinity;
22 for (let i = 0; i < n; i++) {
23 const r = (i > 0 ? Math.abs(e[i - 1]) : 0) + (i < n - 1 ? Math.abs(e[i]) : 0);
24 lo = Math.min(lo, d[i] - r);
25 hi = Math.max(hi, d[i] + r);
26 }
27 const tol = 1e-13 * Math.max(Math.abs(lo), Math.abs(hi), 1e-300);
28 const out: number[] = [];
29 for (let k = 0; k < Math.min(count, n); k++) {
30 let a = lo;
31 let b = hi;
32 while (b - a > tol) {
33 const mid = 0.5 * (a + b);
34 if (countBelow(d, e, mid) > k) b = mid;
35 else a = mid;
36 if (mid === a && mid === b) break;
37 }
38 out.push(0.5 * (a + b));
39 }
40 return out;
41}