back to table · edit · history · where entries came from · files · download
6531 bytes, as of the version from 2026-10-25 13:49 (current). Recorded here, not run.
"""Eigenvalues of the clamped plate -- numberdb.org/T488.
Run with SageMath (or passagemath-symbolics) and numberdb installed:
sage -pip install numberdb mpmath # once
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # fill the existing draft only
sage -python checks.py # independent checks, offline
sage -python checks.py --stored # check the database values
Set NUMBERDB_API_KEY for authenticated access to the unpublished draft.
--fill sends entries, and does not offer, review or publish the table.
The first 30 eigenvalues on the radius-one disk, counted with multiplicity,
are the fourth powers of cross-product roots below 51/4. This cutoff falls
between complete eigenspaces. All roots are certified with Arb signs at exact
rational endpoints. mpmath supplies guesses only. See certificates.md for the
root-count argument and checks.py for the independent radial ODE computation.
The 180-bit guard leaves more than 110 decimal digits in every returned ball.
"""
import os
import sys
from functools import lru_cache
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
import mpmath as mp
WORKING_GUARD = 180
CUTOFF = QQ(51)/4
MAX_ORDER = 9
def rational_decimal(text):
f = Fraction(str(text))
return QQ(f.numerator)/f.denominator
def sign(x):
assert x.is_finite(), 'nonfinite sign enclosure'
if x > 0:
return 1
if x < 0:
return -1
raise ArithmeticError('sign not certified')
def bessel(C, kind, m, x):
z = getattr(C(x), 'bessel_' + kind)(m)
assert z.real().is_finite() and z.imag().is_finite() and z.imag().contains_zero()
return z.real()
def cross(C, m, x):
return (bessel(C, 'J', m, x)*bessel(C, 'I', m+1, x)
+ bessel(C, 'J', m+1, x)*bessel(C, 'I', m, x))
def bisect_certified(f, left, right, width):
sl, sr = sign(f(left)), sign(f(right))
assert sl != sr
while right-left > width:
mid = (left+right)/2
sm = sign(f(mid))
if sm == sl:
left = mid
else:
right = mid
return left, right
def j_brackets(C, m):
# There are no J_m zeros in (0,1]. On [1,infinity), sqrt(x) J_m(x)
# has zero spacing > 2 for every integer m>=0 (Sturm comparison).
# Hence a unit mesh and certified signs count *all* the zeros to 18.
brackets = []
left = QQ(1)
previous = sign(bessel(C, 'J', m, left))
assert previous == 1
for stop in range(2, 19):
right = QQ(stop)
current = sign(bessel(C, 'J', m, right))
if current != previous:
brackets.append(bisect_certified(
lambda x: bessel(C, 'J', m, x), left, right, QQ(1)/10**14))
left, previous = right, current
assert brackets and brackets[-1][0] > CUTOFF
return brackets
def approximate_root(m, left, right, working):
with mp.workdps(working):
a = mp.mpf(str(left.numerator()))/mp.mpf(str(left.denominator()))
b = mp.mpf(str(right.numerator()))/mp.mpf(str(right.denominator()))
def f(x):
return mp.besselj(m+1,x)+mp.besselj(m,x)*mp.besseli(m+1,x)/mp.besseli(m,x)
root = mp.findroot(f, (a+(b-a)/5, a+2*(b-a)/5))
assert a < root < b
return mp.nstr(root, working-10)
@lru_cache(maxsize=None)
def spectrum(digits=100):
C = ComplexBallField(numberdb.bits(digits, losing=WORKING_GUARD))
R = RealBallField(C.precision())
modes = []
evidence = []
eps = QQ(1)/10**(digits+12)
for m in range(MAX_ORDER+1):
js = j_brackets(C, m)
for k in range(1, len(js)):
# A root is between the kth and (k+1)st J_m zeros. If the kth
# zero already exceeds the cutoff, all later roots do too.
if js[k-1][0] > CUTOFF and k > 1:
break
a, b = js[k-1][1], js[k][0]
q = rational_decimal(approximate_root(m, a, b, digits+40))
lo, hi = q-eps, q+eps
assert a < lo < hi < b
assert sign(cross(C,m,lo)) != sign(cross(C,m,hi))
root = R((lo+hi)/2) + R(0,(hi-lo)/2)
derivative = (2*bessel(C,'J',m,root)*bessel(C,'I',m,root)
- cross(C,m,root)/root)
sign(derivative)
assert root.is_finite() and root.rad() < abs(root.center())/2**350
evidence.append((m,k,lo,hi,js[k-1],js[k]))
if lo > CUTOFF:
break
assert hi < CUTOFF
eigenvalue = root**4
modes.append((eigenvalue,m,k))
if m == MAX_ORDER:
assert evidence[-1][0:2] == (MAX_ORDER,1)
assert evidence[-1][2] > CUTOFF
# gamma_{m,k} increases with m; m=9's first root excludes all m>=9.
modes.sort(key=lambda row: row[0].center())
for a,b in zip(modes,modes[1:]):
assert a[0] < b[0], 'distinct mode enclosures overlap'
rows=[]
for value,m,k in modes:
rows.extend([(value,m,k)]*(1 if m==0 else 2))
assert len(rows)==30
return rows,evidence
class ClampedPlate(numberdb.Generator):
table = 'T488'
files = ('generate.py', 'checks.py', 'certificates.md')
parameters = ('domain','n')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for n in range(1,31):
yield {'domain':'unit-disk','n':n}
def value(self,params,digits):
assert params['domain']=='unit-disk'
n=int(params['n'])
if not 1<=n<=30:
raise ValueError('This generator covers the first 30 modes.')
value,m,k=spectrum(digits)[0][n-1]
result={'number':value}
if n==1:
result['comment']='The fundamental mode is axisymmetric.'
if n>1 and spectrum(digits)[0][n-2][1:]==(m,k):
result['equals']='HREF{Eigenvalues_of_the_clamped_plate#unit-disk,%d}'%(n-1)
return result
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator=ClampedPlate()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
print(generator.publish(message='Fill 30 certified clamped-disk eigenvalues'))
else:
report=generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)