Riemann sum not converging

Question

I am trying to find the volume of a hemisphere by numerical integration. I have a set of points equally spaced over x-y plane. I am trying to calculate its volume by trying to evaluate Riemann sum. $$V= \sum_i (\Delta x_i)(\Delta y_i)z_i$$

$$V= \sum_i (\Delta a_i)z_i$$

$$V= \frac{\pi R^2}{N} \sum_i z_i $$

where N is the total number of points.

This is how the points on the sphere look

As I evaluate this sum with increasing number of points, the error should go to zero. But I am observing that it is not at all reducing beyond ~21%. How the error is (not) becoming arbitrarily small

I have used the following python code:

R=10 

volume_approx=[]

error=[]

volume_actual=2/3*np.pi*(R**3)

pts_range=range(10,10000)

for i in pts_range:

    pts = create_uniform_hemisphere(i, R) #(this creates a uniform hemisphere with i points and radius=R)

    base_area=np.pi*(R**2)

    volume=(base_area/i)*np.sum(pts[:,2])

    volume_approx.append(volume)

    err=(volume_actual-volume)/volume_actual

    error.append(100*err)

Can anyone tell me what mistake I am making? Similar behavior is repeated when I try to find volume of a cone, or use randomly spaced points. But this behavior is NOT repeated when I try to calculate volume of a cylinder or a paraboloid with similar code.

Edit: I have used the following code to create the hemisphere:

def create_uniform_hemisphere(n, radius):
    step=2*radius/(math.floor(math.sqrt(n)))
    X,Y = np.mgrid[-radius:radius:step, -radius:radius:step] #create mesh on x-y plane

    X=X.flatten()
    Y=Y.flatten()

    X=np.reshape(X,(X.shape[0],1))
    Y=np.reshape(Y,(Y.shape[0],1))
    pts=np.append(X,Y,axis=1)

    pts=pts[(X.flatten()**2 + Y.flatten()**2)<radius**2] #filter points out of circle

    Z=(radius**2-(pts[:,0]**2 + pts[:,1]**2))**0.5 #project onto hemisphere
    Z=np.reshape(Z,(Z.shape[0],1))

    pts=np.append(pts,Z, axis=1)

    return pts

Answer 1

From what you posted (images and code) it seems you are doing roughly the right thing. The big unknown is the function "create_uniform_hemisphere(i, R)".

You are assuming in your code and the formula (by simply adding the $z$-values of your hemisphere and scaling the sum with $\frac{base\_area}i$) that each point of your sample is representing the same area ($\frac{base\_area}i$) in the projection to the $xy$-plane.

My guess is that this is not true!

It's easy to put 9, 16 (or any square number of) points into a square and have each point represend the same area of the base. It's much harder to do that with 17 points. Now try to do that for a circle!

My suggestion: Make a 2D plot that plots the points (x=pts[0],y=pts[1], z=0) in the $xy$-plane, similiar to the 3d plot that you posted. Look at this plot directly from above ($z$-direction). Do that for for a few small $i$ (try 20, 25, 30) and see if there are actually 20/25/30 points in the circle (and not just that many points in the square that encompasses the circle).

Next, do it for a big $i$, like 100,000. Try to confirm visually that the points are roughly equally dense in each area. It suspect that they may not be, maybe more dense in the center or something similar.

ADDED: I can't really understand the added code at the moment, but the most obvious potential problem case is that I assume that the length of the returned array pts is smaller than the $i$ you use in the main problem. I assume that from the way you create the array. Please add a check in your main loop to compare $i$ and the lenght of the array pts!

If I'm correct, this would also account for the value of your difference ($21\%)$.

Because if you have exactly $i$ points in the bounding square, you have $\frac{\pi}4i$ points in the circle. $1-\frac{\pi}4=0,2146\ldots.$