back to table · edit · history · where entries came from · files · download
13335 bytes, as of the version from 2026-09-27 14:04 (current). Recorded here, not run.
"""Reproduce the 100-digit upper-branch Falkner-Skan tables.
Install dependencies with:
sage -pip install numberdb python-flint mpmath
Check every stored entry (NUMBERDB_API_KEY is needed to read a draft):
sage -python generate.py
Publish a checked precision refinement:
sage -python generate.py --publish
Compute all three tables offline, without reading or writing the website:
sage -python generate.py --offline results
Taylor shooting is run at two working precisions, orders, steps and cutoffs.
The thermal integral includes its Gaussian asymptotic tail. Midpoints are
propagated, so neither Taylor truncation nor the finite-domain boundary
condition is rigorously enclosed. Rigour is heuristic (agreement-checked).
The initial guesses select the existing upper branch, not the final values.
No parameter ranges or mathematical definitions are changed.
"""
import argparse
import json
import time
from decimal import Decimal, localcontext
from fractions import Fraction
from functools import lru_cache
from types import SimpleNamespace
import numberdb
from pathlib import Path
from flint import arb, arb_poly, arb_series, ctx
from mpmath import mp
def rational(text):
value = Fraction(text)
return arb(value.numerator) / value.denominator
def text(value):
return value.mid().str(ctx.dps - 12, radius=False)
def coefficients(state, beta, order):
flow, velocity, acceleration = state
flow_coefficients = [flow, velocity, acceleration / 2]
velocity_coefficients = [velocity, acceleration]
acceleration_coefficients = [acceleration]
for degree in range(order - 2):
product = sum((flow_coefficients[index] * acceleration_coefficients[degree - index]
for index in range(degree + 1)), arb(0))
square = sum((velocity_coefficients[index] * velocity_coefficients[degree - index]
for index in range(degree + 1)), arb(0))
next_acceleration = (beta * square - product - (beta if degree == 0 else 0)) / (degree + 1)
acceleration_coefficients.append(next_acceleration)
velocity_coefficients.append(next_acceleration / (degree + 2))
flow_coefficients.append(next_acceleration / ((degree + 2) * (degree + 3)))
return arb_poly(flow_coefficients)
def integrate(shear, options, moments=False):
beta = rational(options.beta)
step = rational(options.step)
count = int(Fraction(options.limit) / Fraction(options.step))
if count * Fraction(options.step) != Fraction(options.limit):
raise ValueError('Domain length must be a multiple of the step')
state = [arb(0), arb(0), arb(shear)]
displacement, momentum, integral_flow = arb(0), arb(0), arb(0)
prandtls = options.prandtls if moments else []
thermal = {prandtl: arb(0) for prandtl in prandtls}
minimum_acceleration = arb(shear)
for index in range(count):
polynomial = coefficients(state, beta, options.order)
velocity = polynomial.derivative()
acceleration = velocity.derivative()
if moments:
displacement += (1 - velocity).integral()(step)
momentum += (velocity * (1 - velocity)).integral()(step)
integral_series = arb_series(polynomial.integral().coeffs(), prec=options.order)
for prandtl in prandtls:
factor = rational(prandtl)
integrand = (-factor * integral_series).exp().integral()
thermal[prandtl] += (-factor * integral_flow).exp() * arb_poly(integrand.coeffs())(step)
thermal[prandtl] = thermal[prandtl].mid()
integral_flow = (integral_flow + polynomial.integral()(step)).mid()
displacement, momentum = displacement.mid(), momentum.mid()
state = [polynomial(step).mid(), velocity(step).mid(), acceleration(step).mid()]
if not all(value.is_finite() for value in state) or abs(state[1]) > 10:
raise ArithmeticError('Shooting trajectory left the upper-branch neighborhood')
minimum_acceleration = min(minimum_acceleration, state[2])
if not moments:
return state[1] - 1
heat = {}
for prandtl in prandtls:
factor = rational(prandtl)
tail = (-factor * integral_flow + factor * state[0] ** 2 / 2).exp()
tail *= (arb.pi() / (2 * factor)).sqrt() * (state[0] * (factor / 2).sqrt()).erfc()
heat[prandtl] = {'hartree': text(1 / (thermal[prandtl] + tail)),
'wedge': text(1 / ((thermal[prandtl] + tail) * (2 - beta).sqrt())),
'gaussian_tail': text(tail)}
return {'shear_hartree': text(arb(shear)), 'shear_wedge': text(arb(shear) / (2 - beta).sqrt()),
'momentum_hartree': text(momentum), 'momentum_wedge': text(momentum * (2 - beta).sqrt()),
'displacement_hartree': text(displacement), 'heat': heat,
'momentum_identity_residual': text(arb(shear) - beta * displacement - (1 + beta) * momentum),
'endpoint_velocity_residual': text(state[1] - 1),
'endpoint_acceleration': text(state[2]), 'minimum_sampled_acceleration': text(minimum_acceleration)}
DEFAULT_TABLE = 'T417'
SEEDS = {
'-1/6': '0.172091702711', '-1/8': '0.271743686685', '0': '0.469599988361',
'1/4': '0.731940848513', '1/3': '0.802125592789', '1/2': '0.927680039837',
'2/3': '1.03890348316', '3/4': '1.09044156217', '4/5': '1.12026765738',
'1': '1.23258765682',
}
THERMAL_BETAS = ('0', '-1/6', '1/3', '1/2', '2/3', '1')
PRANDTLS = ('1/8', '1/6', '1/4', '1/3', '1/2', '2/3', '1', '3/2', '2', '3', '4', '6', '8')
SETTINGS = ((160, 140, '1/8', '32'), (210, 180, '1/10', '36'))
def solve(beta, setting):
dps, order, step, limit = setting
ctx.dps, ctx.cap, mp.dps = dps, order, dps
options = SimpleNamespace(beta=beta, order=order, step=step, limit=limit,
prandtls=PRANDTLS if beta in THERMAL_BETAS else ())
guess = mp.mpf(SEEDS[beta])
def residual(shear):
return mp.mpf(text(integrate(mp.nstr(shear, dps), options)))
started = time.monotonic()
root = mp.findroot(residual, (guess, guess + mp.mpf('1e-10')),
tol=mp.mpf(10) ** (-(dps - 35)), maxsteps=14)
result = integrate(mp.nstr(root, dps), options, moments=True)
result.update(beta=beta, settings=list(setting), elapsed_seconds=time.monotonic() - started)
return result
def quantities(result):
values = {key: result[key] for key in ('shear_hartree', 'shear_wedge',
'momentum_hartree', 'momentum_wedge')}
for prandtl, row in result['heat'].items():
for normalisation in ('hartree', 'wedge'):
values['heat:' + prandtl + ':' + normalisation] = row[normalisation]
return values
def check_pair(first, second):
if first['beta'] != second['beta'] or quantities(first).keys() != quantities(second).keys():
raise ValueError('Profiles do not describe the same parameters')
with localcontext() as context:
context.prec = 240
differences = {}
for field, value in quantities(first).items():
other = Decimal(quantities(second)[field])
value = Decimal(value)
if not value.is_finite() or not other.is_finite() or min(value, other) <= 0:
raise ArithmeticError('Nonpositive or nonfinite quantity: ' + field)
relative = abs(value / other - 1)
if relative >= Decimal('1e-110'):
raise ArithmeticError('Insufficient convergence: ' + field + ': ' + str(relative))
differences[field] = str(relative)
for result in (first, second):
for field in ('momentum_identity_residual', 'endpoint_velocity_residual',
'endpoint_acceleration'):
value = Decimal(result[field])
if not value.is_finite() or abs(value) >= Decimal('1e-110'):
raise ArithmeticError('Failed boundary or momentum control: ' + field)
if Decimal(result['minimum_sampled_acceleration']) < Decimal('-1e-110'):
raise ArithmeticError('Not the monotone upper branch')
if result['beta'] == '0':
for normalisation in ('hartree', 'wedge'):
relative = abs(Decimal(result['heat']['1'][normalisation]) /
Decimal(result['shear_' + normalisation]) - 1)
if relative >= Decimal('1e-110'):
raise ArithmeticError('Reynolds analogy failed')
return differences
@lru_cache(maxsize=10)
def profile(beta):
first, second = (solve(beta, setting) for setting in SETTINGS)
check_pair(first, second)
return second
class FalknerSkan(numberdb.Generator):
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
def __init__(self, table=DEFAULT_TABLE):
if table not in ('T417', 'T418', 'T419'):
raise ValueError('Choose T417, T418 or T419')
self.table = table
self.parameters = (('beta', 'branch', 'prandtl', 'normalisation') if table == 'T417'
else ('beta', 'branch', 'normalisation'))
def enumerate(self):
betas = THERMAL_BETAS if self.table == 'T417' else tuple(SEEDS)
for beta in betas:
for prandtl in PRANDTLS if self.table == 'T417' else (None,):
for normalisation in ('hartree', 'wedge'):
params = dict(beta=beta, branch='upper', normalisation=normalisation)
if prandtl is not None:
params['prandtl'] = prandtl
yield params
def from_profile(self, params, result, digits=100):
if not 1 <= digits <= 100:
raise ValueError('This generator has been checked only up to 100 significant digits')
if params['branch'] != 'upper' or params['beta'] not in SEEDS:
raise ValueError('Unsupported branch or pressure gradient')
if params['normalisation'] not in ('hartree', 'wedge') or result['beta'] != params['beta']:
raise ValueError('Profile parameters do not match the requested entry')
if self.table == 'T417':
value = result['heat'][params['prandtl']][params['normalisation']]
else:
field = 'shear' if self.table == 'T418' else 'momentum'
value = result[field + '_' + params['normalisation']]
with localcontext() as context:
context.prec = 240
return format(Decimal(value), '.' + str(digits - 1) + 'e')
def value(self, params, digits):
if params['branch'] != 'upper' or params['beta'] not in SEEDS:
raise ValueError('Unsupported branch or pressure gradient')
return self.from_profile(params, profile(params['beta']), digits)
def offline(directory):
directory.mkdir(parents=True, exist_ok=True)
profiles, diagnostics = {}, {}
source_hash = __import__('hashlib').sha256(Path(__file__).read_bytes()).hexdigest()
manifest_path = directory / 'source-sha256.txt'
if manifest_path.exists() and manifest_path.read_text().strip() != source_hash:
raise ValueError('Refuse to reuse cached results from another source')
manifest_path.write_text(source_hash + '\n')
for beta in SEEDS:
pair = []
for index, setting in enumerate(SETTINGS):
path = directory / (beta.replace('/', '_') + '-' + str(index) + '.json')
if path.exists():
result = json.loads(path.read_text())
if result['beta'] != beta or result['settings'] != list(setting):
raise ValueError('Cached settings changed')
else:
print('Computing', beta, setting, flush=True)
result = solve(beta, setting)
path.write_text(json.dumps(result, indent=2) + '\n')
pair.append(result)
diagnostics[beta] = check_pair(*pair)
profiles[beta] = pair[1]
print('Checked', beta, max(Decimal(value) for value in diagnostics[beta].values()), flush=True)
for table in ('T417', 'T418', 'T419'):
generator = FalknerSkan(table)
entries = [{'params': params, 'number': generator.from_profile(params, profiles[params['beta']])}
for params in generator.enumerate()]
(directory / (table + '-entries.json')).write_text(json.dumps(entries, indent=2) + '\n')
(directory / 'diagnostics.json').write_text(json.dumps(diagnostics, indent=2) + '\n')
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--table', choices=('T417', 'T418', 'T419'), default=DEFAULT_TABLE)
mode = parser.add_mutually_exclusive_group()
mode.add_argument('--offline', type=Path)
mode.add_argument('--publish', action='store_true')
options = parser.parse_args()
if options.offline:
offline(options.offline)
return
generator = FalknerSkan(options.table)
if options.publish:
print(generator.publish(message='100-digit agreement-checked Falkner-Skan refinement'))
else:
report = generator.verify(sample=None)
print(report)
if not report.ok:
raise SystemExit(1)
if __name__ == '__main__':
main()