back to table · edit · history · where entries came from · files · download
1751 bytes, as of the version from 2026-10-30 04:30 (current). Recorded here, not run.
# Python 3 with SciPy: a double-precision example for rational 0 <= n < 5.
# Eight significant digits are displayed; no rigorous error bound is claimed.
# Tighten tolerances and check convergence, especially near n = 5.
from fractions import Fraction
from math import isfinite, pi
from scipy.integrate import solve_ivp
def pressure(index, rtol=2e-13, atol=2e-15, start=1e-4, limit=1e5):
n_exact = Fraction(index)
if not 0 <= n_exact < 5:
raise ValueError('Require 0 <= n < 5.')
n = float(n_exact)
if not 0 <= n < 5:
raise ValueError('Index too close to an endpoint for double precision.')
if not 0 < start < limit:
raise ValueError('Require 0 < start < limit.')
# Regular central series through xi^4 and its derivative.
theta = 1 - start**2 / 6 + n * start**4 / 120
slope = -start / 3 + n * start**3 / 30
def rhs(xi, state):
y, v = state
# A zero-density extension only locates the surface event;
# fractional powers are evaluated only on nonnegative arguments.
density = 1.0 if n == 0 else max(y, 0.0)**n
return [v, -2 * v / xi - density]
def surface(xi, state):
return state[0]
surface.terminal = True
surface.direction = -1
sol = solve_ivp(rhs, (start, limit), [theta, slope], method='DOP853',
events=surface, rtol=rtol, atol=atol)
if sol.status != 1:
raise RuntimeError('Surface not reached: ' + sol.message)
v_surface = sol.y_events[0][0, 1]
result = 1 / (4 * pi * (n + 1) * v_surface**2)
if not isfinite(result) or result <= 0:
raise ArithmeticError('Invalid pressure constant.')
return result
if __name__ == '__main__':
print(f"{pressure('13/4'):.8g}")