容易看出$f(1/2)=s(1/2)=1/2$，因此$(1/2,1/2)$是曲线与圆的一个交点。令
$$g(x)=\frac{1}{2}-\sqrt{\frac{1}{16}-(x-\frac{1}{4})^2},$$
则$g(x)$表示圆在直线$y=1/2$以下的部分（包括圆与该直线的交点）。从图中看，另一个交点在$(0,1/2)$内，且在$g(x)$上，将其横坐标记作a，再令
$$h(x)=f(x)-g(x).$$
则需要计算$\int_a^{1/2}h(x)dx$。这可以分解为三个子问题：找出$a$，计算$\int_a^{1/2}f(x)dx$，以及计算$\int_a^{1/2}g(x)dx$。
由于f处处不可导，也不是Lipschitz函数，牛顿法、不动点法等常见的方法都不可用，不过f是连续的，因此仍可以用二分法来计算a。
根据Wiki的信息（可以参考比如这个页面），$\int_a^{1/2}f(x)dx$可以用微积分基本定理求解，f的原函数满足
$$I(x)=\int_{0}^x f(x)dx=\begin{cases} \frac{I(2x)}{4}+\frac{x^2}{2} & 0\le x\le \frac{1}{2} \\ \frac{1}{2}-I(1-x) & \frac{1}{2}\le x\le 1\end{cases}$$
可以递归处理。
最后，$\int_a^{1/2}g(x)dx$不难解析计算，不过为了减少计算量，可以记两交点分别为A、B，圆心为O，先算出AB与x轴所围成的梯形的面积，再减去圆上AB劣弧所形成的弓形的面积。
而弓形的面积等于扇形OAB的面积减去三角形OAB的面积。最终结果是$0.11316017$。
注：本题所涉及的积分知识是我国《高等数学》课程的标准内容。
注2：以下代码为Python 3，作浮点数除法比Python 2更直观一些。
import math
epsilon=0.0000000001
bound=100
f=lambda x: sum(abs(x*(1 < n)-round(x*(1 < n)))/(1 < n) for n in range(bound))
g=lambda x: 0.5-math.sqrt(0.5*x-x*x)
def I(x):
if x < epsilon:
return 0
elif x <= 0.5:
return I(2*x)/4+x*x/2
else:
return 0.5-I(1-x)
a,b=0,0.5
for i in range(bound):
m=(a+b)*0.5
if f(m)>g(m):
b=m
else:
a=m
b=0.5
ya,yb=g(a),g(b)
r=1/4
d=math.sqrt((a-b)*(a-b)+(ya-yb)*(ya-yb))
theta=math.acos(1-0.5*d*d/(r*r))
s=(ya+yb)*(b-a)*0.5-r*r*(theta-math.sin(theta))*0.5
print("%.8f" %(I(b)-I(a)-s))
|
