mpmath--mpmath
wehub-resource-sync
64 行
1.9 KiB
Python
64 行
1.9 KiB
Python
'''
|
||
This script calculates the constant in Gerver's solution to the moving sofa
|
||
problem.
|
||
|
||
See Finch, S. R. "Moving Sofa Constant." §8.12 in Mathematical Constants.
|
||
Cambridge, England: Cambridge University Press, pp. 519-523, 2003.
|
||
'''
|
||
|
||
from mpmath import cos, sin, pi, quad, findroot, mp
|
||
|
||
mp.prec = 113
|
||
|
||
eqs = [lambda A, B, φ, θ: (A*(cos(θ) - cos(φ)) - 2*B*sin(φ)
|
||
+ (θ - φ - 1)*cos(θ) - sin(θ) + cos(φ) + sin(φ)),
|
||
lambda A, B, φ, θ: (A*(3*sin(θ) + sin(φ)) - 2*B*cos(φ)
|
||
+ 3*(θ - φ - 1)*sin(θ) + 3*cos(θ) - sin(φ) + cos(φ)),
|
||
lambda A, B, φ, θ: A*cos(φ) - (sin(φ) + 0.5 - 0.5*cos(φ) + B*sin(φ)),
|
||
lambda A, B, φ, θ: ((A + pi/2 - φ - θ) - (B - (θ - φ)*(1 + A)/2
|
||
- 0.25*(θ - φ)**2))]
|
||
A, B, φ, θ = findroot(eqs, (0, 0, 0, 0))
|
||
|
||
def r(α):
|
||
if 0 <= α < φ:
|
||
return 0.5
|
||
if φ <= α < θ:
|
||
return (1 + A + α - φ)/2
|
||
if θ <= α < pi/2 - θ:
|
||
return A + α - φ
|
||
return B - (pi/2 - α - φ)*(1 + A)/2 - (pi/2 - α - φ)**2/4
|
||
|
||
s = lambda α: 1 - r(α)
|
||
|
||
def u(α):
|
||
if φ <= α < θ:
|
||
return B - (α - φ)*(1 + A)/2 - (α - φ)**2/4
|
||
return A + pi/2 - φ - α
|
||
|
||
def du(α):
|
||
if φ <= α < θ:
|
||
return -(1 + A)/2 - (α - φ)/2
|
||
return -1
|
||
|
||
def y(α, f):
|
||
if α > pi/2 - θ:
|
||
i = [0, φ, θ, pi/2 - θ, α]
|
||
elif α > θ:
|
||
i = [0, φ, θ, α]
|
||
elif α > φ:
|
||
i = [0, φ, α]
|
||
else:
|
||
i = i = [0, α]
|
||
return 1 - quad(lambda x: f(x)*sin(x), i)
|
||
|
||
y1 = lambda α: y(α, r)
|
||
y2 = lambda α: y(α, s)
|
||
y3 = lambda α: y2(α) - u(α)*sin(α)
|
||
|
||
S1 = quad(lambda x: y1(x)*r(x)*cos(x), [0, φ, θ, pi/2 - θ, pi/2 - φ])
|
||
S2 = quad(lambda x: y2(x)*s(x)*cos(x), [0, φ, θ])
|
||
S3 = quad(lambda x: y3(x)*(u(x)*sin(x) - du(x)*cos(x) - s(x)*cos(x)),
|
||
[φ, θ, pi/4])
|
||
|
||
print(2*(S1 + S2 + S3))
|