一个球被两个垂直的平面切割，分成四份，已知三份体积，求剩下的那份体积？

nyy

https://bbs.emath.ac.cn/forum.ph ... 3&fromuid=14149

这个问题的推广提问

nyy

先简单一点，只是被两个垂直的面切，然后看情况，是否被三个垂直的面切

TSC999

一个圆球，被两个互相垂直的平面切成体积不同的四部分。求最小那一块的体积。有以下公式：




例题：

1. Clear["Global`*"];
   
2. R = 10; Subscript[h, 1] = 3; Subscript[h, 2] = 5;
   
3. V = 2/3 Subscript[h, 1] Subscript[h, 2] Sqrt[R^2 - \!\(\*SubsuperscriptBox[\(h\), \(1\), \(2\)]\) - \!\(\*SubsuperscriptBox[\(h\), \(2\), \(2\)]\)] + 2/3 R^3 ArcTan[(R Sqrt[R^2 - \!\(\*SubsuperscriptBox[\(h\), \(1\), \(2\)]\) -         
   
4. \!\(\*SubsuperscriptBox[\(h\), \(2\), \(2\)]\)])/(Subscript[h, 1] Subscript[h, 2])] - 1/3 Subscript[h, 2] (3 R^2 - \!\(\*SubsuperscriptBox[\(h\), \(2\), \(2\)]\)) ArcTan[Sqrt[R^2 - \!\(\*SubsuperscriptBox[\(h\), \(1\), \(2\)]\) -        
   
5. \!\(\*SubsuperscriptBox[\(h\), \(2\), \(2\)]\)]/Subscript[h, 1]] - (1/3) Subscript[h, 1] (3 R^2 - \!\(\*SubsuperscriptBox[\(h\), \(1\), \(2\)]\)) ArcTan[Sqrt[R^2 - \!\(\*SubsuperscriptBox[\(h\), \(1\), \(2\)]\) - \!\(\*SubsuperscriptBox[\(h\), \(2\), \(2\)]\)]/Subscript[h, 2]]
   
6. N[V]

复制代码

运行结果：

\(10 \sqrt{66}-\frac{1375}{3} \arctan\left(\sqrt{\frac{22}{3}}\right)+\frac{2000}{3} \arctan\left(2\sqrt{\frac{22}{3}}\right)-291 \arctan\left(\frac{\sqrt{66}}{5}\right)\)

152.346

我做这个题，大概是20年前的事情了。记得在本论坛上贴出过，帖子名称记不得了。

yuange1975

球和圆的差距不大。

球的15.268
圆的15.385

yuange1975

前天爬山手机上写代码调试出来，今天有空整理了一下。

懒得写分析文章了，代码很简单，有兴趣的自己看代码吧。小小的打一个广告，有兴趣的可以加入我的知识星球《袁哥的技术天地》，主要是网络安全、数学、投资等方面的一些经验的随意分享。



两互相垂直平面切一球，切成四部分，依次体积为abcd。已知abc,求d=?

这题情形现实中还是很常见的，切球形西瓜、土豆等。

就是图1的圆的问题推广成球的问题。球半径R,两垂直平面交点P坐标(R*x,R*y)。

n为精度位数，可以设置成需要的求解精度要求。
求得如图1的数据d=15.2684403090431785。


abcd按逆时针顺序排列，就是分别是坐标体系里的第一二三四象限。

有三个方程：
3(a+b)/(pi*R^3)=y^3-3y+2
3(b+c)/(pi*R^3)=-x^3+3x+2

r=sqrt(1-x*x-y*y)
rx=sqrt(1-x*x)
ry=sqrt(1-y*y)

a=f(x,y)=(2*x*y*r+2*acos(x/rx*y/ry)+(x**3-3*x)*acos(y/rx)+(y**3-3*y)*acos(x/ry)/3*R^3

求d的方程
3(a+b+c+d)=4*pi*R^3

P点坐标(x*R,y*R)(x^2+y^2<=1)切的第一块a块的体积公式f(x,y)。b=f(-x,y)和c=f(-x,-y)、d=f(x,-y)都有这样的公式，但是这个公式太复杂。a+b、b+c可以消去一些复杂项，公式变得非常简单。就取一个复杂的a就行，另外两个用简单的方程。

解这3个方程组得到xyR，带入就得到d=4/3*pi*R^3-(a+b+c)。



Python 代码如下：

1. # 两互相直平面切一球，切成四部分，依次体积为abcd。
   
2. # 已知abc,求d=?
   
3. # copy by yuange 2023.9.20
   
4. # n 精度位数
   
5. 
   
6. from  mpmath import *
   
7. 
   
8. n=20
   
9. 
   
10. mp.prec=128
    
11. mp.dps=n+1
    
12. rn="\r\n"
    
13. 
    
14. 
    
15. '''
    
16. 
    
17. 3(a+b)/(piR^3)=Y^3-3Y+2
    
18. 3(b+c)/(piR^3)=(-X^3+3X+2)
    
19. Y^3-3Y+2=(a+b)/(b+c)*(-X^3+3X+2)
    
20. a=R^3*vv(x,y)
    
21. 
    
22. '''
    
23. 
    
24. stepn=16
    
25. num=1000
    
26. 
    
27. mpf1=mpf("1")
    
28. mpf0=mpf("0")
    
29. er=pow(mpf("0.1"),n)
    
30.      
    
31. def vv(x,y):
    
32.     r=sqrt(1-x*x-y*y)
    
33.     v0=(2*x*y*r+2*acos(x*y/sqrt(r*r+x*x*y*y))+(x**3-3*x)*acos(y/sqrt(r*r+y*y))+(y**3-3*y)*acos(x/sqrt(r*r+x*x)))/3
    
34.     return v0
    
35. 
    
36. def qiu(a,b,c):
    
37.     step=mpf1/stepn
    
38.     y=0-mpf1
    
39.     stop=0
    
40.     sign=0
    
41.     if c>=b :
    
42.        y=mpf0
    
43.        stop=1
    
44.     fx=lambda x :-x**3+3*x+2-mpf1*(b+c)/(a+b)*(y**3-3*y+2)
    
45.     for i in range(num):
    
46.         fy=(y**3-3*y+2)*(b+c)/(a+b)-2
    
47.         if abs(fy)>2 :
    
48.            y+=step
    
49.            continue
    
50.         x=findroot(fx,0)
    
51.         #print(rn,"x,y=",x,y)        
    
52.         if x*x+y*y > 1+1e-5 :
    
53.            y+=step
    
54.            if y > 1 :break
    
55.            continue
    
56.         a1=vv(x,y)
    
57.         R1=3*(a+b)/pi/(y**3-3*y+2)
    
58.         R2=a/a1      
    
59.         #print(rn,"R1,R2=",R1,R2,R1-R2,step,a1)        
    
60.         if abs(R1-R2) < er :
    
61.            R=cbrt(R1)
    
62.            d=4*pi*R1/3-a-b-c
    
63.            #print(rn,"ok!X,Y,R,d=",x*R,y*R,R,d)
    
64.            return x,y,R,d
    
65.            break
    
66.         if (R1-R2)>0 :
    
67.            y-=step           
    
68.            step=step/stepn
    
69.            #print(rn,"step=",step)
    
70.            continue      
    
71.         y+=step         
    
72.         if y>stop+0.1 :
    
73.            break
    
74. 
    
75. x,y,r,d=qiu(12,20,25)
    
76. print(rn,"x,y,r,d=",(x,y,r,d))
    
77. x,y,r,d=qiu(20,25,d)
    
78. print(rn,"x,y,r,d=",(x,y,r,d))

复制代码

yuange1975

最初v(x,y)是使用的数值积分，不能直接解方程，所以等于是自己写了一个自己求解方程的代码。
得到了v(x,y)的解析表达式后就可以直接用find root 解方程组了，代码就非常简单了。

部分方程得不到解析解，可以数值积分等的可以使用那个框架自己解方程。

Python代码和结果如下，代码和结果刚好一屏。

1. # 两互相垂直平面切一球，切成四部分，依次体积为abcd。
   
2. # 已知abc,求d=?
   
3. # copy by yuange 2023.9.20
   
4. # n 精度位数
   
5. from  mpmath import *
   
6. n=20
   
7. mp.dps=n+1
   
8. 
   
9.      
   
10. def v(x,y):
    
11.     r=sqrt(1-x*x-y*y)
    
12.     rx=sqrt(1-x*x)
    
13.     ry=sqrt(1-y*y)
    
14.     v0=(2*x*y*r+2*acos(x/rx*y/ry)+(x**3-3*x)*acos(y/rx)+(y**3-3*y)*acos(x/ry))/3
    
15.     return v0
    
16. def qiu(a,b,c):
    
17.     f1=lambda x,y:y**3-3*y+2+(x**3-3*x-2)*(a+b)/(b+c)
    
18.     f2=lambda x,y:a*pi/v(x,y)-3*(a+b)/(y**3-3*y+2)
    
19.     x,y=findroot([f1,f2],[0,0])
    
20.     d=4*(a+b)/(y**3-3*y+2)-(a+b+c)
    
21.     R=cbrt((a+b+c+d)/pi*3/4)   
    
22.     return (x,y,R,d)   
    
23. print(qiu(12,20,25))

复制代码

yuange1975

这个公式非常不错，不过形式还是有点不太好。

球半径R。过球心，做垂直两垂直面的垂面，切出来横截面，取两垂直面分别水平和垂直方向，交点P的坐标(x*R,y*R)。限制P位于圆上或者圆内，有x^2+y^2<=1


切下来第一象限那块体积：

V(x,y)=1/3*(2*x*y*r+2*acos(x/rx*y/ry)+(x**3-3*x)*acos(y/rx)+(y**3-3*y)*acos(x/ry))*R^3

其中
r=sqrt(1-x*x-y*y)
rx=sqrt(1-x*x)
ry=sqrt(1-y*y)

这样的公式形式更好更漂亮，也更适用。

V(x,y)对于x或者y的负值也适用，并且第二三四象限的块体积分别为，V(-x,y)、V(-x,-y)、V(x,-y)。

这样更通用，可以一次直接适应各种大小面积分布。

得到公式很重要，还有公式的形式也很重要。
一个是公式怎么更通用，比如这个里面的xy范围不要局限于都大于0。
一个是公式各自变量怎么更独立，比如这个xy取坐标对于R的比例，而不是直接坐标，可以把R独立出来。
还有怎么写成更简洁的形式，比如这个公式里把r、rx、ry独立表示出来，原公式整体就会显得简洁很多。

DyannaMenphina

TSC999 发表于 2023-8-20 09:59
一个圆球，被两个互相垂直的平面切成体积不同的四部分。求最小那一块的体积。有以下公式：

您好，能请教您一下关于这个积分的思路吗？

DyannaMenphina

TSC999 发表于 2023-8-20 09:59
一个圆球，被两个互相垂直的平面切成体积不同的四部分。求最小那一块的体积。有以下公式：

您好，我想请教一下在计算这个问题的时候，那个包含了(x-sin(x))的那一项积分的思路。