The background to this problem is that I have a computational code that is used to model fluid flow on a Cartesian computational grid. I have a disk of material that is in rotation about a central point and the material is gravitationally bound.

My problem is that I had been calculating the gravitational force at the center of each grid cell which will have an error as the force of gravity is inversely proportional to r^2.

I'm trying to make my calculations more accurate by taking a volume integral of the gravitational force as the errors can be fairly large when r is fairly small.

The force due to gravity is given by Fgrav=GM/(x^2+y^2), where G and M are constants.

The fact that I'm working on a cartesian grid (which I have no choice about) makes the double integral of Fgrav with x and y as the variables hard and I'm not sure how to go about this one, the limits of the integral will the the vertices of a square basically. Has anyone any suggestions on how to tackle it?

Thanks in advance for any advice.