4

I want to calculate the volume (V) of a part of sphere, which is the result of intersetion of the sphere with three palnes (x=0, y=0 and z=1.5). I am using R-Language and this is my code. I tried 2 different methods using cartesian and polar coordinates. Both of them deliver negative answers. 

##  Compute the Volume between 3 planes x=0, y=0 and z=1.5 and a sphere
library("pracma", lib.loc="~/R/win-library/3.1")
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1.5 ) # here the function(x,y)  is subtracted with -1.5 which represents the plane z=1.5 
xmin <- 0
xmax <- 2
ymin <- 0 
ymax <- function(x) (sqrt(4 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # Integral over sector is not so exact

# Exact Volume from AutoCAD V=0.3600



##  Volume of the sphere: use polar coordinates
f0 <- function(x, y) (sqrt(4 - x^2 - y^2)-1.5)  # for x^2 + y^2 <= 4 the f(x,y) means z changes between zmin=1 and zmax= sqrt(4-x^2-y^2)
fp <- function(t, r) r * f0(r*cos(t), r*sin(t))
quad2d(fp, 0, pi/2, 0, 2, n = 101)  # -0.523597

The correct answer is V= 0.3600 . Can anyone give a me hint, please?

Cheers

Mohamad Reza Salehi Sadaghiani
  • 295
  • 1
  • 10

2 Answers2

5

Your X-Y region of integration covers areas where f(x,y)-1.5 is negative, as well as positive. The intersection of your sphere with the line z=1.5 is a circle of radius sqrt(7/4) (using Pythagoras), so adjusting your limits appropriately, you get:

library(pracma)
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1.5 ) # here the function(x,y)  is subtracted with -1.5 which represents the plane z=1.5 
xmin <- 0
xmax <- sqrt(7/4)
ymin <- 0 
ymax <- function(x) (sqrt(7/4 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # Integral over sector is not so exact
# [1] 0.3599741

Pretty close to what you're expecting.

Miff
  • 7,486
  • 20
  • 20
  • That was very helpfull. But when I want to calculate the volume between 3 planes x=1,y=1,z=1 and same sphere with radius 2. I changed the limits (x from 1 to sqrt(3) and y from 1 to sqrt(3-x^2)). It gives wrong result (V= 0.0330) . The right one is V= 0.0152. – Mohamad Reza Salehi Sadaghiani Jun 12 '15 at 14:28
  • No, it's not. At least if I do understand your integration domain. The correct numerical value would be 0.26113536... ! Correct the upper y limit to `sqrt(4-x^2)` -- as it should be -- and *integral2* will return 0.251140 . You may try to improve the accuracy by setting the `reltol` parameter, but I think you will get better results by using polar coordinates. – Hans W. Jun 13 '15 at 13:39
  • I have used AutoCAD and the exact Volume of a Piece between 3 Planes x=1,y=1 and z=1 and a sphere x^2+y^2+z^2=4 is V=0.0152. – Mohamad Reza Salehi Sadaghiani Jun 15 '15 at 10:39
0

I have solved it for the sphere with r=2 and the volume of intersection between sphere shell and 3 planes x=1,y=1 and z=1.

##  Compute the Volume between 3 planes x=1.0, y=1.0 and z=1.0 and a sphere
library(pracma)
f <- function(x, y)  (sqrt(4 -x^2 - y^2) - 1 ) # here the function(x,y) is  subtracted with -1.5 which represents the plane z=1.5 
xmin <- 1
xmax <- sqrt(2)
ymin <- 1 
ymax <- function(x) (sqrt(3 - x^2))
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q   # 
# [1] 0.01520549
# Exact Volume from AutoCAD: 0.0152
Mohamad Reza Salehi Sadaghiani
  • 295
  • 1
  • 10