back to table · edit · history · where entries came from · files · download
8116 bytes, as of the version from 2026-10-25 12:01 (current). Recorded here, not run.
"""Neumann eigenvalues of the classical planar domains -- numberdb.org/T485.
Run with SageMath (or passagemath with scipy and mpmath):
sage -pip install numberdb scipy mpmath
sage -python generate.py # verify all stored entries
sage -python generate.py --fill # fill this draft, key in environment
sage -python checks.py # independent local checks
The first 20 eigenvalues include the constant mode and all multiplicities.
Integer quadratic forms enumerate the integrable polygons. scipy's double
precision Bessel derivative roots select low disk and sector modes; every
selected root is recomputed with mpmath at 130 and 160 decimal digits, then
bracketed by exact rational endpoints using Arb evaluations of J'_m.
The endpoint bracket, rather than the mpmath point estimate, is squared.
The precision guard is 192 bits. The exact root bracket has radius
10^(-digits-12), leaving over ten guard digits after squaring. This is measured
by checks.py. Index identification is agreement checked, not interval certified;
the declared overall rigour is therefore heuristic (agreement-checked).
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
import mpmath as mp
from scipy.special import jnp_zeros
TABLE = 'T485'
COUNT = 20
WORKING_GUARD = 192
ORDER_BOUND = 40
RADIAL_BOUND = 20
DOMAINS = (
'unit-square', 'rectangle-1x2', 'unit-disk', 'half-disk',
'quarter-disk', 'sector-60', 'equilateral-triangle',
'right-isosceles-triangle', 'hemiequilateral-triangle',
)
CIRCULAR = {'unit-disk': 1, 'half-disk': 1, 'quarter-disk': 2, 'sector-60': 3}
def rational_decimal(text):
"""Convert decimal text to an exact rational without float coercions."""
from fractions import Fraction
x = Fraction(str(text))
return QQ(x.numerator) / QQ(x.denominator)
def derivative(order, x, C):
z = C(x)
if order == 0:
return -z.bessel_J(1).real()
return (z.bessel_J(order - 1) - z.bessel_J(order + 1)).real() / 2
@lru_cache(None)
def root_enclosure(order, k, digits):
# mpmath counts the zero at the origin for order zero; scipy does not.
index = k + (order == 0)
with mp.workdps(digits + 30):
low_precision = mp.besseljzero(order, index, derivative=1)
with mp.workdps(digits + 60):
high_precision = mp.besseljzero(order, index, derivative=1)
assert abs(low_precision - high_precision) < mp.mpf(10)**(-digits - 20)
center = rational_decimal(mp.nstr(high_precision, digits + 20))
radius = QQ(1) / (ZZ(10)**(digits + 12))
a, b = center - radius, center + radius
bits = numberdb.bits(digits, losing=WORKING_GUARD)
R, C, I = RealBallField(bits), ComplexBallField(bits), RealIntervalField(bits)
fa, fb = derivative(order, R(a), C), derivative(order, R(b), C)
assert fa.is_finite() and fb.is_finite()
assert (fa > 0 and fb < 0) or (fa < 0 and fb > 0), (order, k)
root = R(I(a, b))
assert root.is_finite() and root > 0 and root.accuracy() > numberdb.bits(digits + 9, losing=0)
return root
@lru_cache(None)
def layouts():
"""Mode labels; integer/rational sorting for all polygonal domains."""
result = {}
for domain in DOMAINS:
if domain in CIRCULAR:
continue
rows = []
for j in range(COUNT + 1):
for k in range(COUNT + 1):
if domain in ('right-isosceles-triangle', 'hemiequilateral-triangle') and k > j:
continue
if domain in ('equilateral-triangle', 'hemiequilateral-triangle'):
q = QQ(16) * (j*j + j*k + k*k) / 9
elif domain == 'rectangle-1x2':
q = QQ(j*j) + QQ(k*k) / 4
else:
q = QQ(j*j + k*k)
rows.append((q, j, k))
rows.sort()
# Any omitted j or k is >= COUNT+1; even the smallest coefficient
# (1/4 in the rectangle) then exceeds the last retained coefficient.
assert rows[COUNT - 1][0] < QQ((COUNT + 1)**2) / 4
result[domain] = [('pi', q, j, k) for q, j, k in rows[:COUNT]]
roots = {m: jnp_zeros(m, RADIAL_BOUND) for m in range(ORDER_BOUND + 1)}
for domain, step in CIRCULAR.items():
candidates = [(0.0, 0, 0, 0)]
for m in range(0, ORDER_BOUND + 1, step):
multiplicity = 2 if domain == 'unit-disk' and m > 0 else 1
for k, root in enumerate(roots[m], 1):
candidates.extend((float(root)**2, m, k, copy) for copy in range(multiplicity))
candidates.sort()
cutoff = candidates[COUNT - 1][0]
# DLMF 10.21.3: the first positive derivative root for m>0 is >m.
assert cutoff < (ORDER_BOUND + 1)**2
assert all(float(roots[m][-1])**2 > cutoff + 1 for m in range(0, ORDER_BOUND + 1, step))
result[domain] = [('zero',)] + [('bessel', m, k) for _, m, k, _ in candidates[1:COUNT]]
return result
class NeumannEigenvalues(numberdb.Generator):
table = TABLE
parameters = ('domain', 'n')
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
files = ('generate.py', 'checks.py')
def enumerate(self):
for domain in DOMAINS:
for n in range(1, COUNT + 1):
yield {'domain': domain, 'n': n}
def number(self, params, digits):
mode = layouts()[params['domain']][int(params['n']) - 1]
if mode[0] == 'zero' or (mode[0] == 'pi' and mode[1] == 0):
return ZZ(0)
if mode[0] == 'pi':
R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
return R(mode[1]) * R.pi()**2
return root_enclosure(mode[1], mode[2], digits)**2
def value(self, params, digits):
result = {'number': self.number(params, digits)}
domain, n = params['domain'], int(params['n'])
mode = layouts()[domain][n - 1]
# These exact mode identities give the corresponding Dirichlet
# fundamental eigenvalue. They are also checked against stored T390.
fundamental = {
'unit-square': QQ(2), 'rectangle-1x2': QQ(5)/4,
'equilateral-triangle': QQ(16)/3,
'right-isosceles-triangle': QQ(5),
'hemiequilateral-triangle': QQ(112)/9,
}
if mode[0] == 'pi' and mode[1] == fundamental.get(domain):
result['equals'] = 'HREF{Dirichlet_eigenvalues_of_the_classical_planar_domains#%s,1}' % domain
elif mode == ('bessel', 0, 1):
result['equals'] = 'HREF{Dirichlet_eigenvalues_of_the_classical_planar_domains#half-disk,1}'
else:
# Repeat multiplicities and shared exact modal values without
# mistaking numerical closeness for equality.
def identity(item):
if item[0] == 'zero' or (item[0] == 'pi' and item[1] == 0):
return ('zero',)
return item[:2] if item[0] == 'pi' else item
wanted = identity(mode)
for other in DOMAINS:
for i, candidate in enumerate(layouts()[other], 1):
if (other, i) == (domain, n):
return result
if identity(candidate) == wanted:
result['equals'] = 'HREF{#%s,%s}' % (other, i)
return result
return result
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
generator = NeumannEigenvalues()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill the checked first 20 Neumann eigenvalues of nine domains.'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)