back to table · edit · history · where entries came from · files · download
5495 bytes, as of the version from 2026-10-25 13:49 (current). Recorded here, not run.
"""Independent checks for Eigenvalues of the clamped plate, numberdb.org/T488.
Run: sage -python checks.py [--stored]
The stored mode needs NUMBERDB_API_KEY while T488 is a draft.
No special-function evaluations enter the independent ODE determinant.
"""
import os
import sys
import json
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import mpmath as mp
from generate import ClampedPlate, spectrum
def rational(s):
q=Fraction(str(s));return QQ(q.numerator)/q.denominator
def written_interval(s):
d=Decimal(s)
unit=QQ(10)**int(d.as_tuple().exponent)
q=rational(s)
return q-unit,q+unit
def ode_basis(lam,m,parity,R):
# Regular solutions of D_m^2 u=lam*u: sum c_j*r^(m+2j).
# c_j = lam*c_(j-2)/(16*j*(j-1)*(m+j)*(m+j-1)).
term=R(1);v=R(1);der=R(m+2*parity)
last=parity
for j in range(parity+2,181,2):
term=term*lam/(16*j*(j-1)*(m+j)*(m+j-1))
v+=term;der+=(m+2*j)*term;last=j
j=last+2
nextterm=term*lam/(16*j*(j-1)*(m+j)*(m+j-1))
ratio=lam/(16*(j+2)*(j+1)*(m+j+2)*(m+j+1))
assert ratio>0 and ratio<QQ(1)/2
assert ratio*QQ(m+2*(j+2))/(m+2*j)<QQ(1)/2
# Both subsequent ratios decrease, so a geometric series bounds tails.
v+=R(0,2*nextterm.upper())
der+=R(0,2*(m+2*j)*nextterm.upper())
return v,der
def ode_det(lam,m,R):
if lam==0:
# The recurrence has zero tails, avoiding its strict positive test.
return R(2)
a,ad=ode_basis(lam,m,0,R);b,bd=ode_basis(lam,m,1,R)
return a*bd-b*ad
def ode_det_mp(lam,m):
cols=[]
for parity in (0,1):
term=mp.mpf(1);v=mp.mpf(1);d=mp.mpf(m+2*parity)
for j in range(parity+2,1000,2):
term*=lam/(16*j*(j-1)*(m+j)*(m+j-1))
v+=term;d+=(m+2*j)*term
if j>30 and abs(term)*(m+2*j)<mp.eps:
break
else:raise ArithmeticError('series did not converge')
cols.append((v,d))
return cols[0][0]*cols[1][1]-cols[0][1]*cols[1][0]
def strict_sign(x):
assert x.is_finite()
assert x.rad()<abs(x.center())/2
assert x>0 or x<0
return 1 if x>0 else -1
def read_rows():
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
doc=numberdb.table('T488');numbers=doc['Numbers']['unit-disk']
return {int(n):e['number'] for n,e in numbers.items() if isinstance(e,dict) and 'number' in e}
g=ClampedPlate()
return {p['n']:to_text(g.value(p,100)['number'],digits=100) for p in g.enumerate()}
def main():
R=RealBallField(768)
# Nonzero known controls exercise the actual series evaluator at lambda=0.
for m in range(10):
assert ode_det_mp(mp.mpf(0),m)==2
for parity in (0,1):
# For lambda=0 the exact solution is r^(m+2*parity).
assert (1*(m+2)-m*1)==2
# Positive-lambda control: known Taylor coefficient of determinant at 0.
# det=2-lambda/(8*(m+1)*(m+2)*(m+3))+O(lambda^2).
for m in range(10):
with mp.workdps(60):
got=mp.diff(lambda z:ode_det_mp(z,m),0)
expected=-mp.mpf(1)/(8*(m+1)*(m+2)*(m+3))
assert mp.almosteq(got,expected)
print('Controls: ODE boundary determinant at zero = 2; first derivative correct for orders 0..9.')
rows=read_rows();assert sorted(rows)==list(range(1,31))
modes,certificates=spectrum()
checked={};min_accuracy=10000;low_last=None
for n in range(1,31):
text=rows[n];lo,hi=written_interval(text);ball,m,k=modes[n-1]
assert len(Decimal(text).as_tuple().digits)==100
assert lo>0 and hi-lo<QQ(1)/10**94
assert R(lo)<ball<R(hi)
min_accuracy=min(min_accuracy,int(ball.accuracy()))
if (m,k) not in checked:
left=ode_det(R(lo),m,R);right=ode_det(R(hi),m,R)
assert strict_sign(left)!=strict_sign(right)
roots=[]
for prec in (140,180):
with mp.workdps(prec):
center=mp.mpf(text)
root=mp.findroot(lambda z:ode_det_mp(z,m),(center*mp.mpf('0.999'),center*mp.mpf('1.001')))
rootq=rational(mp.nstr(root,prec-10))
assert lo<rootq<hi
roots.append(rootq)
assert abs(roots[0]-roots[1]) < abs(roots[1])*QQ(1)/10**125
checked[m,k]=text
else:assert checked[m,k]==text
if low_last is not None and modes[n-2][1:]!=(m,k):
assert low_last < lo
low_last=hi
print('PASS n=%d m=%d k=%d %s'%(n,m,k,text),flush=True)
# Published fundamental eigenvalue 104.363, a lower precision external check.
assert abs(rational(rows[1])-rational('104.363'))<rational('0.0005')
print('Independent ODE signs and 140/180-digit roots agree for all 30 entries (%d distinct modes).'%len(checked))
print('Minimum returned-ball accuracy (bits):',min_accuracy)
print('Characteristic root certificates:',len(certificates))
print('Serialized values:',json.dumps(rows,sort_keys=True))
print('Root certificates:',json.dumps([{'m':m,'k':k,'lo':str(lo),'hi':str(hi),'j_left':[str(x) for x in a],'j_right':[str(x) for x in b]} for m,k,lo,hi,a,b in certificates]))
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':numberdb.configure(api_key=sys.stdin.read().strip())
main()