Numerical integration over a disk with polar coordinates

[uzi kiko]

In my job, I was given the task of calculating a force that operates an ultrasound transmitter on a receiver. The calculation is made by assuming that each point on the transmitter is a small transmitter and integration should be made on the surface of the transmitter.
Since the transmitter is disc-shaped, it would seem natural to use polar coordinates. But the results I obtained were not relatively accurate compared to the analytical results.
I enclose here a simple Peyton code that compares numeric integration on a unit-sized disk (the area is pi), by a summation of rings, in polar coordinates, and in Cartesian coordinates.
The script results:
rings= [3.1415926535897114, 2000]
polar= [3.143163449916506, 4002000]
cartzian=[3.141529000085212, 4000000]
pi= 3.141592653589793

import numpy as np
import math

# rings
dr = 1.0/2000.0
dang = 2*np.pi/2000.0
radius = 1s=0
r=dr
i=0
while (r<radius): 
    ang=0
    ring = 0
    s += np.pi * (2*r*dr -dr*dr)
    r+=dr
    i+=1
print("rings= " + str([s, i]))# polar
s=0
r=dr
i=0
while (r<radius): 
    ang = 0
    dang_ = dang
    while (ang<2*np.pi):
        i+=1
        ang+=dang_ 
    #ang -= (np.random.rand() - 0.5) * dang * 2;
    factor = 2*np.pi/ang
    ang *= factor
    s += ang * r * dr
    r+=dr
print("polar= " + str([s, i]))

# cartesian
dx = 1.0/1000.0
x=-1
y=-1
s1 = 0
i=0
while (x<1):
    y= -1   
    while (y<1):       
        if math.sqrt(x**2+y**2)<=radius:
            s1+=dx*dx
        y+=dx
        i+=1       
    x+=dx   
        
print("cartzian=" + str([s1,i]))
print("pi= " + str(np.pi))

It can be easily seen that in the calculation in polar coordinates the error is already in the third digit after the dot. To reach the same level of accuracy of the Cartesian coordinates, I need to reduce the resolution to almost 1000 times (which makes the calculation impractical in a commercial system).
Does anyone recognize where my mistake is in calculating the area in polar coordinates?

 

[pbuk]

When you are integrating over a surface you want to choose the size and shape of your elements $\delta A$ so that

1. the function is nearly constant over each $\delta A$ individually (in order to be accurate), and also so that
2. the difference between the function and the approximation is similar over all $\delta A$s (otherwise you waste time generating a high accuracy in some places which is lost in others). You also need to consider the implication of your choice of coordinates on the effects of roundoff error.

Unless your objective is simply to calculate the area of a shape, you won't get much idea of how a method performs in general by using a constant function: you need to use a function that is more like the ones you need to support.

I haven't looked at your methods in detail, but at a quick glance:

• rings is 'cheating' because it is using the fact that the function is constant over a ring. This is fine if it is true for your real world functions of course (i.e. they have infinite rotational symmetry), otherwise you would also need to calculate around the ring. This is similar to a polar calculation with the order of the coordinates swapped.
• I think polar is losing accuracy because of point 2 above - you are calculating tiny areas near the centre compared to much larger ones near the boundary. Because your elements are not square there is an error in each one.
• cartesian works fine: for your simple function the approximation is exact except on elements that straddle the border, and so there is little loss of accuracy.

 

[uzi kiko]

Thanks a lot!
Exactly what I was looking for.

I enclose here the code I changed according to your answer so that now the calculation of the area using a polar integral is correct up to the fifth digit after the dot, and about half the steps for the same calculation in Cartesian coordinates.

s=0
r=dr
i=0
diff = 0
da = np.pi*dr*(2.0*radius-dr)/2000.0
while (r<radius): 
    ang = 0
    N = 2000*(2*r-dr)/(2*radius-dr)
    N+=diff
    i0 = 0
    while (i0<N):
        i+=1
        s+=da
        i0+=1
    diff = N-i0
    r+=dr
print("polar= " + str([s, i]))

 

[pasmith]

In my job, I was given the task of calculating a force that operates an ultrasound transmitter on a receiver. The calculation is made by assuming that each point on the transmitter is a small transmitter and integration should be made on the surface of the transmitter.

Is there some reason you cannot use a library function from https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html? This has the advantage that the issues raised by @pbuk will already have been considered, and the algorithm extensively tested and debugged.

 

[uzi kiko]

Thanks, @pasmith.
Our system is based on c++. I used here Python code just to illustrate the problem.
BTW:
I am not sure that I agree with @pbuk regarding his point about the polar errors. Although there is a difference in resolution close to the center relative to the resolution at the end of the circle. The size of the area is also different.

 

[pasmith]

Thanks, @pasmith.
Our system is based on c++. I used here Python code just to illustrate the problem.

In that case I would suggest investigating the GNU Scientific Library.

If you want to write your own routines, then Numerical Recipes is your starting point.