Here is a simplied version of
Werner's code, assumed
diameter = 1
I like dimensionless version because you can get a sense of results.
Example, with diameter of 1, sphere volume = Pi/6 ≈ 0.523598775598 ≈ ½
You can sense how big the hole is ...
Also, A, B is limited to inside a unit circle.
Code:
00 { 83-Byte Prgm }
01▸LBL "RB1"
02 MVAR "A"
03 MVAR "B"
04 MVAR "V"
05 1
06 RCL "A"
07 X↑2
08 -
09 RCL "B"
10 X↑2
11 -
12 SQRT
13 LSTO "K"
14 RCL "A"
15 RCL× "B"
16 STO× ST Y
17 RCL÷ "K"
18 ATAN
19 -
20 3
21 ÷
22 RCL "A"
23 RCL "B"
24 XEQ 10
25 RCL- "V"
26 RCL "B"
27 RCL "A"
28▸LBL 10
29 RCL÷ "K"
30 ATAN
31 X<>Y
32 STO× ST Y
33 X↑2
34 3
35 -
36 6
37 ÷
38 ×
39 -
40 END
Example, to get sphere - hole, diameter = 20, hole = 2x2
Solve with A=0.1, B=0.1, we have hole = 0.00996658845699, about 2% of sphere
Pi/6 - V = 0.513632187141
Scale up by D^3 => sphere - hole = 0.513632187141 * 20^3 = 4109.05749713
Update:
Perhaps we should compare volume against maximum possible volume
Solve with sphere - 4 caps, or simply solve with A = B = √(0.5), we get
V(max) ≈ 0.402001836518
So, above hole ≈ 2½ % of maximum possible volume.
(06-11-2020 06:53 PM)Albert Chan Wrote: [ -> ]Solve with A=0.1, B=0.1, we have hole = 0.00996658845699, about 2% of sphere
For small square hole, we can estimate with simple formula
For unit diameter (d=1), with square hole, b=a
XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> hv2(b = 0.1) // just to confirm square hole formula
→ 0.00996658845699
XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)
Since taylor series all have even powers, let dimensionless x = b^2
XCas> pade(x - x^2/3 - 7*x^3/90, x, 3,3)
→ (17*x^2-30*x)/(7*x-30)
XCas> r(x) := x*(17*x-30)/(7*x-30)
XCas> r(0.1^2) // rough estimate for small 0.1×0.1 hole
→ 0.009966588
70698
(06-16-2020 12:35 PM)Albert Chan Wrote: [ -> ]For unit diameter (d=1), with square hole, b=a
XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)
hv2 ≈ b^2 * √(1 - (2b^2)/3) = b^2 * √(1 - (1-k^2)/3)
Remove b=a requirement: hv2 ≈ a*b * √(1 - (a^2+b^2)/3)
Remove d=1 requirement:
hv2 ≈ a*b * √(d^2 - (a^2+b^2)/3)
It is not as accurate hv2 pade approximation, but likely good enough
Previous example, d = 1, a = b = 0.1
lua> d, a, b = 1, 0.1, 0.1
lua> a*b * sqrt(d*d - (a*a+b*b)/3)
0.009966610925150703
Hole volume with 5 good digits, over-estimated.
Last term (estimated average height) is equivalent to
RMS(d, d, sqrt(d*d-a*a-b*b))
Repalce RMS with arithemetic mean, we get the under-estimated hole volume.
lua> a*b * (2*d + sqrt(d*d-a*a-b*b)) / 3
0.00996649831220389
This is not as good as RMS version, does not reduce computation, thus not recommended.
To be clear, my son's integral solves for the volume left behind in the sphere after the hole has been removed. I have listed three volumes below. I wasn't sure if I was clear on this.
d = 12
a = 2
b = 2
Volume of sphere solid = 904.779
Volume of sphere with hole removed = 857.0638
Volume of hole removed = 47.715
Thoughts?
(11-05-2023 04:28 PM)DM48 Wrote: [ -> ]d = 12
a = 2
b = 2
...
Volume of sphere solid = 904.779
Volume of sphere with hole removed = 857.0638
Volume of hole removed = 47.715
sphere volume = pi*d^3/6, looks good --> we check removed hole only.
This is a good test of my recent estimate, lo < (average height) < hi
lua> d, a, b = 12, 2, 2
lua> lo = (2*d + sqrt(d*d-a*a-b*b))/3
lua> hi = sqrt(d*d - (a*a+b*b)/3)
lua> a*b*lo, a*b*hi
47.5492050529208 47.55347866700536
Your hole volume is outside this range.
sphere (with hole) volume formula is incorrect.
Perhaps you can show steps to derive the formula?
Here is Derive6 calculated hole volume
#1: r := 6
#2: a := 2
#3: b := 2
#4: 8·∫(∫(√(r·r - x·x - y·y), x, 0, a/2), y, 0, b/2)
Shift-Enter (or click [✔≈]) → 47.55262983
Comment:
8 = 2^3, we could scale-up all 3 dimensions, to remove factor 8
#5 d := 2*r
#6: ∫(∫(√(d·d - x·x - y·y), x, 0, a), y, 0, b) → 47.55262983
This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.
(11-05-2023 08:01 PM)DM48 Wrote: [ -> ]This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.
I'd say you can bring down the size of that image to about 10-15% of it's current size, it's silly the way it is now.
(11-05-2023 08:27 PM)rprosperi Wrote: [ -> ] (11-05-2023 08:01 PM)DM48 Wrote: [ -> ]This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.
I'd say you can bring down the size of that image to about 10-15% of it's current size, it's silly the way it is now.
I have tried to resize it with img=##x## but it does not display.
(11-05-2023 09:52 PM)DM48 Wrote: [ -> ]I have tried to resize it with img=##x## but it does not display.
Not sure what's going on, as Joe said maybe there are limits. I wasn't suggesting you make it unusable, it just looked a bit silly the way it was. Feel free to restore it if that's the only way it's useful to see the details.
(11-05-2023 04:16 PM)Albert Chan Wrote: [ -> ] (06-16-2020 12:35 PM)Albert Chan Wrote: [ -> ]For unit diameter (d=1), with square hole, b=a
XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)
hv2 ≈ b^2 * √(1 - (2b^2)/3) = b^2 * √(1 - (1-k^2)/3)
Remove b=a requirement: hv2 ≈ a*b * √(1 - (a^2+b^2)/3)
Remove d=1 requirement: hv2 ≈ a*b * √(d^2 - (a^2+b^2)/3)
To show hole volume is over-estimated, we check sign of errors
lua> taylor(hv2 - b^2 * √(1 - (2b^2)/3), b, 10)
- b^6/45 - 23*b^8/945 - 143*b^10/5670 + b^12*order_size(b)
Negative under-estimated error → RMS formula over-estimated hole volume.
Quote:Last term (estimated average height) is equivalent to RMS(d, d, sqrt(d*d-a*a-b*b))
Repalce RMS with arithemetic mean, we get the under-estimated hole volume.
height of rect cap (both of them), h = 1-k // 1 is diameter, k is corner height
mean([1,1,k]) = 1 - h/3 = (1-h) + 2/3*h
MEAN formula under-estimated hole volume, if rect cap mean height > 2/3*h
XCas> ratio := (hv2/(b*b)-k) / (1-k)
XCas> taylor(ratio, b, 8, polynom)
2/3 + 4/45*b^2 + 5/63*b^4 + 23/252*b^6 + 1997/16632*b^8
Small hole, h and b^2 about the same size: k = 1-h = sqrt(1-2*b^2) ≈ 1-b^2
Ratio converge faster (smaller coefficients), if we use variable h.
XCas> taylor(ratio(b=sqrt((1-(1-h)^2)/2)), h, 6, polynom)
2/3 + 4/45*h + 11/315*h^2 + 1/84*h^3 + 25/8316*h^4
Compare mean height ratio with spherical cap.
V = pi*h/6 * (3*a^2 + h^2)
A = pi*a^2
ratio = V/A/h = 1/2 + (h/a)^2/6 > 1/2
Code:
function hv2(d,a,b) -- sphere, rectanglar hole removed volume
local d2, a2, b2 = d*d, a*a, b*b
local k = sqrt(d2 - a2 - b2)
d = a*b*k - d*d2*atan((a*b)/(d*k))
return (2*d + a*(3*d2-a2)*atan(b/k) + b*(3*d2-b2)*atan(a/k)) / 6
end
-- fast estimate: m1 < hv2 < m2
function hv2_m2(d,a,b) return a*b * sqrt(d*d - (a*a+b*b)/3) end
function hv2_m1(d,a,b) return a*b * (2*d + sqrt(d*d-a*a-b*b))/3 end
Although estimation formula assumed square hole, inequality work with rectangular hole too.
For estimate, hv2_m2(d,a,b) is recommended: cheaper to compute, better accuracy.
Link to derivation of formula.
\( \int_{0}^{a} \int_{0}^{\sqrt{d^2-y^2-b^2}}\sqrt{d^2-y^2-x^2}-b\text{ }dx\text{ }dy\text{ }+\int_{a}^{d}\int_{0}^{\sqrt{d^2-y^2}}\sqrt{d^2-y^2-x^2}\text{ }dx\text{ }dy \)
In[842]:= Clear[a, b, q, r, x]
d = 12
a = 2
b = 2
q = NIntegrate[
Integrate[
Sqrt[d^2 - y^2 - x^2] - b, {x, 0, Sqrt[d^2 - y^2 - b^2]}], {y, 0,
a}]
r = NIntegrate[
Integrate[Sqrt[d^2 - y^2 - x^2], {x, 0, Sqrt[d^2 - y^2]}], {y, a, d}]
q + r
Out[848]= 857.2260543948448
a=4
b=2
d=12
Out[849]=811.0424359298954
Hi, DM48
Nice to know there is another way to get volume:
sphere - hole = (x-strip - hole) + (2 spherical caps)
Since we have direct formula for spherical cap, 2nd integral can be skipped
2 caps = 2 * pi/3*h^2*(3r-h) = pi/12*(2h)^2*(3d-2h) = pi/12*(d-a)^2*(2d+a)
10 DEF FNH(X) @ N=N+1 @ FNH=SQRT(H-X*X) @ END DEF
20 DEF FNZ(Y) @ H=D*D-Y*Y @ FNZ=INTEGRAL(0,FNH(B),P,FNH(IVAR)-B) @ END DEF
30 INPUT "D,A,B = ";D,A,B @ N=0 @ P=1E-5
40 I1=INTEGRAL(0,A,P,FNZ(IVAR)) @ I2=PI/12*(D-A)^2*(2*D+A)
50 DISP I1;"+";I2;"=";I1+I2,N @ GOTO 30
>run
D,A,B = 12,2,2
176.547653773 + 680.678408277 = 857.22606205 960
D,A,B = 12,2,10
17.2338339713 + 680.678408277 = 697.912242248 992
D,A,B = 12,10,2
662.307552043 + 35.6047167408 = 697.912268784 1888
Numerically, thinner, less curvy strip (A ≤ B) converge better.
For comparison, below patch directly get (-hole) + sphere
>20 DEF FNX(Y) @ H=D*D-Y*Y @ FNX=INTEGRAL(0,B,P,FNH(IVAR)) @ END DEF
>40 I1=-INTEGRAL(0,A,P,FNX(IVAR)) @ I2=PI/6*D*D*D
>run
D,A,B = 12,2,2
-47.552634005 + 904.778684234 = 857.226050229 225
D,A,B = 12,2,10
-206.866458199 + 904.778684234 = 697.912226035 465
D,A,B = 12,10,2
-206.866453454 + 904.778684234 = 697.91223078 465
Hi, DM48
Double integral reduced to single. Nice!
I simplified a bit, using
Arc SOHCAHTOA method
10 DEF FNZ(Y) @ N=N+1 @ Y=D*D-Y*Y @ FNZ=Y*ACOS(B/SQR(Y)) @ END DEF
20 INPUT "D,A,B = ";D,A,B @ N=0 @ P=1E-5
30 I1=INTEGRAL(0,A,P,FNZ(IVAR))/2 @ C=SQR(D*D-A*A-B*B)
40 I2=PI/12*(D-A)^2*(2*D+A) - B/4*(A*C+(D*D-B*B)*ATAN(A/C))
50 DISP I1;"+";I2;"=";I1+I2,N @ GOTO 20
>run
D,A,B = 12,2,2
200.09879005 + 657.127264322 = 857.2260543
72 15
D,A,B = 12,2,10
82.547121343 + 615.365121066 = 697.912242
409 31
D,A,B = 12,10,2
764.410058495 + -66.4978178412 = 697.91224
0654 31
Below is direct formula for I1 (c is corner height, see code)
>a*(3*d^2-a^2)*pi/12 - a*b*c/12
858.633041958
>res - b*(3*d^2+b^2)*atan(a/c)/12
785.468340817
>res - a*(3*d^2-a^2)*atan(b/c)/6
615.99486321
>res + d^3*atan(a*b/(c*d))/3
764.410059992