back to table · edit · history · where entries came from · files · download
6744 bytes, as of the version from 2026-10-25 12:37 (current). Recorded here, not run.
"""Roots of the frequency equations of the Euler-Bernoulli beam -- numberdb.org/T486.
Run with SageMath (or a compatible passagemath installation):
sage -pip install numberdb
sage -python generate.py # verify all stored entries, read only
sage -python generate.py --fill # write entries to this existing draft
Reading this private draft requires NUMBERDB_API_KEY. The --fill option writes
numbers; it does not publish a table or offer it for review.
The six classical endpoint pairs are those of Khasawneh and Segalman,
arXiv:1811.03222v3, Table 5. We store z=beta*L, not its square or fourth power.
n counts ALL modes with multiplicity. Free-free starts with two exact zeros;
pinned-free starts with one. All other roots are positive and simple.
For positive-root index k, we isolate the cantilever root in
((k-3/4)*pi,(k-1/4)*pi), the clamped-clamped root in
((k+1/4)*pi,(k+3/4)*pi), and the clamped-pinned root in
(k*pi,(k+1/4)*pi). Every initial bracket is checked for opposite signs and
for a derivative interval excluding zero. Rational bisection with Arb signs
then certifies the enclosure. No approximate root is wrapped in a ball.
The generator uses 96 guard bits beyond the 100-digit target. Bisection stops
at absolute width 10^(-digits-12), giving at least twelve extra decimal
places before serialization. For n=1..12 in all six pairs, checks.py
measured a minimum relative-accuracy lower bound of 374 bits (over 112
decimal digits) and independently checked the boundary determinants.
There are no omitted roots between the initial brackets: for the cosine
families, outside these neighborhoods |cos(z)| >= 1/sqrt(2), while sech(z)
is smaller once z >= 3*pi/4. The first cantilever neighborhood contains its
first root; before it cos(z) is positive. For clamped-clamped, cos(z)cosh(z)
decreases from 1 on (0,pi/2) since tan(z)>tanh(z), and cos(z) is negative
from pi/2 to 3*pi/2. For clamped-pinned, tan(z)-tanh(z) is strictly
increasing on every positive branch, is positive on (0,pi/2), and cannot
vanish on a branch where tan(z)<0. It has exactly one root between k*pi
and k*pi+pi/4 for every k>=1.
"""
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.real_mpfi import RealIntervalField
TABLE = 'T486'
BOUNDARIES = ('clamped-free', 'clamped-clamped', 'clamped-pinned',
'pinned-pinned', 'free-free', 'pinned-free')
MAX_N = 12
WORKING_GUARD = 96
def sign(x):
if not x.is_finite():
raise ArithmeticError('A sign enclosure is not finite')
if x > 0:
return 1
if x < 0:
return -1
raise ArithmeticError('The sign is unresolved; increase working precision')
def function(kind, z):
if kind == 'clamped-free':
return z.cos() + 1 / z.cosh()
if kind == 'clamped-clamped':
return z.cos() - 1 / z.cosh()
if kind == 'clamped-pinned':
return z.sin() - z.cos() * z.tanh()
raise ValueError(kind)
def derivative(kind, z):
if kind == 'clamped-free':
return -z.sin() - z.tanh() / z.cosh()
if kind == 'clamped-clamped':
return -z.sin() + z.tanh() / z.cosh()
if kind == 'clamped-pinned':
return z.tanh() * (z.sin() + z.cos() * z.tanh())
raise ValueError(kind)
@lru_cache(maxsize=None)
def positive_root(kind, k, digits):
k = int(k)
if k < 1:
raise ValueError('A positive-root index must be at least one')
precision = numberdb.bits(digits, losing=WORKING_GUARD)
rb = RealBallField(precision)
ri = RealIntervalField(precision)
if kind == 'pinned-pinned':
return k * rb.pi()
if kind == 'clamped-free':
a, b = QQ(k) - QQ(3)/4, QQ(k) - QQ(1)/4
elif kind == 'clamped-clamped':
a, b = QQ(k) + QQ(1)/4, QQ(k) + QQ(3)/4
elif kind == 'clamped-pinned':
a, b = QQ(k), QQ(k) + QQ(1)/4
else:
raise ValueError(kind)
# MPFI endpoints are exact dyadic rationals outside the intended pi bounds.
lo = QQ((a * ri.pi()).lower())
hi = QQ((b * ri.pi()).upper())
sl = sign(function(kind, rb(lo)))
sh = sign(function(kind, rb(hi)))
assert sl * sh == -1, (kind, k, 'no initial sign change')
monotone = sign(derivative(kind, rb(ri(lo, hi))))
assert monotone == sh, (kind, k, 'unexpected derivative sign')
target = QQ(1) / ZZ(10)**(int(digits) + 12)
while hi - lo > target:
mid = (lo + hi) / 2
sm = sign(function(kind, rb(mid)))
if sm == sl:
lo = mid
else:
hi = mid
assert sign(function(kind, rb(lo))) == sl
assert sign(function(kind, rb(hi))) == sh
result = rb(ri(lo, hi))
assert result.is_finite() and QQ(result.rad()) < target
return result
class BeamRoots(numberdb.Generator):
table = TABLE
parameters = ('boundary', 'n')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for boundary in BOUNDARIES:
for n in range(1, MAX_N + 1):
yield {'boundary': boundary, 'n': str(n)}
def value(self, params, digits):
boundary, n = params['boundary'], int(params['n'])
if boundary not in BOUNDARIES or n < 1:
raise ValueError(params)
if boundary == 'free-free':
if n <= 2:
return {'number': ZZ(0), 'comment': (
'Rigid translation, represented by $X(x)=1$.' if n == 1 else
'Rigid rotation, represented by $X(x)=x-L/2$.')}
return {'number': positive_root('clamped-clamped', n-2, digits),
'equals': 'clamped-clamped,%d' % (n-2)}
if boundary == 'pinned-free':
if n == 1:
return {'number': ZZ(0), 'comment':
'Rigid rotation about the pinned end, represented by $X(x)=x$.'}
return {'number': positive_root('clamped-pinned', n-1, digits),
'equals': 'clamped-pinned,%d' % (n-1)}
result = positive_root(boundary, n, digits)
if boundary == 'clamped-free' and n == 1:
return {'number': result, 'comment':
'Fundamental bending mode of the cantilever.'}
return result
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
generator = BeamRoots()
if '--fill' in sys.argv:
print(generator.publish(message='Recompute certified beam roots'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)