"""Steklov eigenvalues of the classical planar domains -- numberdb.org/T489. Run with SageMath (or passagemath) and the numberdb package: sage -pip install numberdb sage -python generate.py # verify every stored entry sage -python generate.py --fill # fill the existing draft, key required sage -python checks.py # independent numerical checks The table holds the first 20 modes for four unscaled domains, including zero and repeating multiplicities. Exact disk modes and the square's bilinear mode are integers. All other values carry Arb enclosures at 100 significant digits. Filling entries does not publish the draft or offer it for review. References: Girouard--Polterovich, doi:10.4171/JST/164, Examples 1.3.1, 4.2.5 and Section 3.1; Girouard--Lagace--Polterovich--Savo, doi:10.1112/S0025579318000414, Lemmas 2.1 and 2.4. Rectangle completeness: a separated eigenfunction is either constant, the bilinear function xy when the half-widths agree, or a trig factor times a hyperbolic factor. On every positive branch, alpha*cot(a*alpha) and -alpha*tan(a*alpha) decrease strictly while alpha*tanh(b*alpha) and alpha*coth(b*alpha) increase strictly. Endpoint limits give one crossing, except the first sine/sinh branch: it has a positive crossing only if b>a. For a=b its zero limit is precisely the bilinear mode, counted separately. Bisection uses rational t=alpha*a/pi; Arb certifies all endpoint signs. For an interval beginning at zero or a pole, rational interior endpoints approach the endpoint until the expected signs are certified. Each retained root's enclosing ball is propagated to sigma. A lower bound for every omitted branch certifies that the first 20 modes are all present. Distinct modes have disjoint balls; the square's two orientations share the same ball. Annulus roots: C_m=(3m/2)*(4^m+1)/(4^m-1); sigma^2-C_m*sigma+m^2/2=0. The smaller root is evaluated by the product relation to avoid cancellation. Each root is double. The radial eigenvalues are 0 and 3/(2*log(2)), simple. Since C_m <= 5m/2, the smaller root is at least m/5. This bounds all omitted angular modes. Rational roots, if any, are detected with integer square roots. Precision: 100 guard bits, with bisection stopped 40 bits before the field precision, retain over 110 relative decimal digits on every nonexact entry. No fixed-precision approximation is ever wrapped and called an enclosure. """ import os import sys from functools import lru_cache from math import isqrt 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 TABLE = 'T489' DOMAINS = ('unit-square', 'rectangle-1x2', 'unit-disk', 'annulus-1-2') LIMIT = 20 GUARD_BITS = 100 def sign(x): assert x.is_finite() if x > 0: return 1 if x < 0: return -1 return 0 def hyperbolic(x, kind): return x.tanh() if kind == 'cosh' else 1 / x.tanh() def root_ball(a, b, trig, hyp, k, field): # Cross-multiplication removes poles, without introducing roots inside # these open positive branches. start = QQ(k) if trig == 'sin' else QQ(k) + QQ(1)/2 stop = start + QQ(1)/2 def f(t): theta = field(t) * field.pi() z = theta * b / a if trig == 'sin': return (theta.cos()-z.tanh()*theta.sin() if hyp == 'cosh' else z.tanh()*theta.cos()-theta.sin()) return (-theta.sin()-z.tanh()*theta.cos() if hyp == 'cosh' else -z.tanh()*theta.sin()-theta.cos()) for power in range(4, 100): epsilon = QQ(1)/ZZ(2)**power lo, hi = start + epsilon, stop - epsilon slo, shi = sign(f(lo)), sign(f(hi)) if slo and shi and slo != shi: break else: raise ArithmeticError(('no sign bracket', a, b, trig, hyp, k)) for _ in range(field.precision()-40): mid = (lo+hi)/2 sm = sign(f(mid)) if not sm: raise ArithmeticError('insufficient precision in root bisection') if sm == slo: lo = mid else: hi = mid assert sign(f(lo)) == slo and sign(f(hi)) == -slo alpha = (field(lo).union(field(hi))) * field.pi()/a return alpha def sort_modes(modes, field): modes.sort(key=lambda row: field(row['number']).mid()) for left, right in zip(modes, modes[1:]): if left['mode'] == right['mode']: continue assert field(left['number']) < field(right['number']), (left,right) return modes def rectangles(a, b, digits): field=RealBallField(numberdb.bits(digits,losing=GUARD_BITS)) modes=[{'number':ZZ(0),'mode':('constant',)}] if a == b: modes.append({'number':1/a,'mode':('bilinear',)}) branch_limit=7 orientations=[(a,b)] if a == b else [(a,b),(b,a)] for ta,hb in orientations: for trig in ('sin','cos'): for hyp in ('cosh','sinh'): for k in range(branch_limit): if trig == 'sin' and hyp == 'sinh' and k == 0 and hb <= ta: continue alpha=root_ball(ta,hb,trig,hyp,k,field) value=alpha*hyperbolic(alpha*hb,hyp) mode=(str(ta),str(hb),trig,hyp,k) row={'number':value,'mode':mode,'alpha':alpha, 'trig':trig,'hyp':hyp,'a':ta,'b':hb,'k':k} modes.extend([row]* (2 if a == b else 1)) sort_modes(modes,field) cutoff=field(modes[LIMIT-1]['number']) for ta,hb in orientations: alpha=field(branch_limit)*field.pi()/ta assert alpha*(alpha*hb).tanh() > cutoff return modes[:LIMIT] def annulus(digits): field=RealBallField(numberdb.bits(digits,losing=GUARD_BITS)) modes=[{'number':ZZ(0),'mode':('constant',)}, {'number':field(3)/(2*field(2).log()),'mode':('radial',)}] max_m=32 for m in range(1,max_m+1): q=ZZ(4)**m c=QQ(3*m)*(q+1)/(2*(q-1)) d=c*c-2*m*m pnum,pden=isqrt(int(d.numerator())),isqrt(int(d.denominator())) if pnum*pnum == d.numerator() and pden*pden == d.denominator(): large=(c+QQ(pnum)/pden)/2 small=QQ(m*m)/(2*large) else: large=(field(c)+field(d).sqrt())/2 small=field(QQ(m*m)/2)/large for branch,value in [('minus',small),('plus',large)]: row={'number':value,'mode':(m,branch),'m':m,'branch':branch} modes.extend([row,row]) sort_modes(modes,field) assert field(QQ(max_m+1)/5)>field(modes[LIMIT-1]['number']) return modes[:LIMIT] @lru_cache(None) def spectrum(domain,digits=100): if domain == 'unit-disk': return [{'number':ZZ(n//2),'mode':('disk',n//2)} for n in range(1,LIMIT+1)] if domain == 'unit-square': return rectangles(QQ(1)/2,QQ(1)/2,digits) if domain == 'rectangle-1x2': return rectangles(QQ(1)/2,QQ(1),digits) if domain == 'annulus-1-2': return annulus(digits) raise ValueError(domain) class SteklovEigenvalues(numberdb.Generator): table=TABLE parameters=('domain','n') type='R' digits=100 rigour='proven' files=('generate.py','checks.py') def enumerate(self): for domain in DOMAINS: for n in range(1,LIMIT+1): yield {'domain':domain,'n':n} def value(self,params,digits): domain=params['domain'];n=int(params['n']) if not 1 <= n <= LIMIT: raise ValueError('This generator enumerates the first 20 modes.') value=spectrum(domain,digits)[n-1]['number'] row={'number':value} if n == 1: row['comment']='The constant eigenfunction gives the simple zero eigenvalue.' row['equals']='HREF{T485#%s,1}'%domain if domain != 'annulus-1-2' else 'unit-square,1' elif domain == 'unit-square' and n == 4: row['comment']=r'The bilinear eigenfunction $(x-1/2)(y-1/2)$ gives the simple eigenvalue $2$.' mode=spectrum(domain,digits)[n-1]['mode'] for earlier in range(1,n): if spectrum(domain,digits)[earlier-1]['mode']==mode: row['equals']='%s,%s'%(domain,earlier) break # The first angular pair of the annulus is an already-tabulated quadratic root. if domain == 'annulus-1-2' and n in (2,3): row['equals']='HREF{Algebraic_numbers_of_degree_2#2,-5,1,1}' if n == 2: row['comment']=r'The first positive eigenvalue is $(5-\sqrt{17})/4$.' if domain == 'annulus-1-2' and spectrum(domain,digits)[n-1]['mode'] == (1,'plus'): row['equals']='HREF{Algebraic_numbers_of_degree_2#2,-5,1,2}' return row if __name__ == '__main__': if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1': os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip() g=SteklovEigenvalues() if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1': print(g.publish(message='Fill 80 Steklov eigenvalues from certified separated modes.')) else: report=g.verify(sample=None) print(report) sys.exit(0 if report.ok else 1)