Numerical integration of gravity at surface of an object

[Ulysees]

Given any distribution of mass, the gravity at any point outside it can be calculated by dividing the mass into infinitesimal masses, and adding up the gravity due to each.

Doing that, I faced a problem on the surface of an object: no matter how small you make the integration step, near the point of interest the distance to the infinitesimal mass gets smaller than the integration step so you get very wrong values. So what is the proper way to integrate gravity for a point at the surface of an object?

Maybe there is a theorem that a hemisphere has this much gravity at its centre, sso you can use the result for infinitesimal masses closer than D from the point of interest on the surface.

Can anyone solve this integral analytically? Total gravity at centre = ?

 

[Ulysees]

 

[merryjman]

I am pretty sure that if you aren't going to do this integration with a computer, the mass distribution will have to exhibit some sort of symmetry, so please be more specific on what your "object" is. For example, Shell Theorem only works if you can approximate your object to be a sphere. I think a hemisphere won't work unless you really look carefully at your indices of integration. With Shell Theorem, you can treat any shell as a point mass located at the shell's center, so the distance from the surface to the shell's center will never be zero.

If you're analyzing a hemisphere, maybe it would be useful to treat the hemisphere as a stack of increasingly larger disks, and find the gravity contribution from each disk, and then integrate over all the disks. That's how Shell works - take a spherical shell, and divide it into "rings," find the gravity force from each ring, and add them all up by integration.

 

[Ulysees]

I am pretty sure that if you aren't going to do this integration with a computer, the mass distribution will have to exhibit some sort of symmetry

The idea is to use numerical integration for an arbitrary distribution that is defined as an equation (eg an ellipsoid like the earth, or an object defined in terms of cubic surface patches). But this fails at the edges between matter and space.

So how do we do the edges? We could increase the resolution locally. That works fine when we're close to the edge (which I've done successfully), but when we're exactly on the edge, it fails.

So we're only stuck when dealing with points on the edge, defined analytically. For these points, we can go like this:

- if current cubic pixel of integration is further than a small radius r, use Newton's law of gravity,

else

- use the analytical formula for a hemisphere of radius r and knowledge of the tangent plane (derived analytically from the equation of the Earth ellipsoid or whatever).

So what is the analytical formula for the gravity at the centre of a hemisphere of radius r, that is the question.

 

[awvvu]

Maybe this could be helpful:

https://www.physicsforums.com/showthread.php?t=225552

I tried doing this problem with a double integral and also a surface integral but I couldn't figure out how to relate dm to the variables of integration.

edit: oh wait, the mass is not uniformly distributed?

 

[Ulysees]

oh wait, the mass is not uniformly distributed?

It is initially.

In future versions of my program it won't be. But if we make r small enough, this will be ok, assuming density does not vary too much locally.

 

[ObsessiveMathsFreak]

Why don't you try some kind of mollified force function.

I take it you are numerically integrating something like
$$\mathbf{F}(\mathbf{x}) = \int \Delta \mathbf{F}(\mathbf{y},\mathbf{x}) d\mathbf{y}$$
Where your "force density" function $$\Delta\mathbf{F}$$ is something like
$$\Delta \mathbf{F}(\mathbf{y},\mathbf{x})=\frac{Gm(\mathbf{x})\rho(\mathbf{y})}{|\mathbf{y}-\mathbf{x}|^2}$$

You could try mollifying the force density to prevent self interaction problems, which are essentially the source of your errors. You could try something like

$$\Delta \mathbf{F}(\mathbf{y},\mathbf{x})=\begin{cases}\frac{4 \pi}{3} G\rho(\mathbf{x})\rho(\mathbf{y})|\mathbf{y}-\mathbf{x}| & |\mathbf{y}-\mathbf{x}|\le k \\ \frac{Gm(\mathbf{x})\rho(\mathbf{y})}{|\mathbf{y}-\mathbf{x}|^2} & |\mathbf{y}-\mathbf{x}|>k. \end{cases}$$
where
$$m(\mathbf{x})=\frac{4 \pi k^3}{3} \rho(\mathbf{x})$$
and k is some characteristic radius which is small, but should be larger than your integration steps.

Or perhaps more simply
$$\Delta \mathbf{F}(\mathbf{y},\mathbf{x})=\begin{cases}0& |\mathbf{y}-\mathbf{x}|\le k \\ \frac{Gm(\mathbf{x})\rho(\mathbf{y})}{|\mathbf{y}-\mathbf{x}|^2} & |\mathbf{y}-\mathbf{x}|>k. \end{cases}$$
Which would avoid self interaction problems entirely, but leave out close by objects.

Let us know if you try these.

 

[Ulysees]

Thanks, I need to think this through when I get some time.

Or perhaps more simply
$$\Delta \mathbf{F}(\mathbf{y},\mathbf{x})=\begin{cases}0& |\mathbf{y}-\mathbf{x}|\le k \\ \frac{Gm(\mathbf{x})\rho(\mathbf{y})}{|\mathbf{y}-\mathbf{x}|^2} & |\mathbf{y}-\mathbf{x}|>k. \end{cases}$$
Which would avoid self interaction problems entirely, but leave out close by objects.

That's close to what I do at present, I can't increase the resolution forever so I pretend there is no mass nearby, which is correct when you're inside matter, but not when you're exactly on the edge, so I end up with a noisy curve for gravity along the surface of the object.

 

[CRGreathouse]

If you're just on the surface, not inside the object, you should be able to find the center of mass by integrating over the volume of the object. This would take three triple integrals in your case, or three double integrals if the density is uniform (since the innermost integrand of the triple integral will be 1, so you can simplify one level).

Once you have that, can't you just calculate gravity as though from a point mass?

 

[Hurkyl]

Once you have that, can't you just calculate gravity as though from a point mass?

Nope. The easiest counterexample to see is to compare the gravitational force from a single point mass, and a pair of point masses with that center of gravity, where one of the points is right next to your test particle.

 

[Hurkyl]

I'm assuming you are treating your original mass as a continuous matter distribution?


You can do some actual analysis and estimate errors. Suppose you are decomposing your mass into little cubes. Then, for each cube, you can compute the smallest and largest densities in that cube. Then, for each component of gravitational force, you should be able to determine the largest and smallest forces that are possible for the cube to exert on your test particle under all possible matter distributions satisfying that density constraint. This way, you have actual hard upper and lower bounds for the value you're calculating.


Another option is to attempt to subtract off the singularity -- before you attempt to do numerical integration, you analyze the integrand and separate it into a singular and a nonsingular part... hopefully where the singular part is easy to do with symbolic integration.

For example, if I wanted to integrate

$$\frac{1}{\sqrt{|\sin x|}}$$

I would recognize that the singluarity at the origin looks like $1 / \sqrt{|x|}$, and I would write the integral as

$$\int_{-1}^1 \frac{1}{\sqrt{|\sin x|}} \, dx = \int_{-1}^1 \frac{1}{\sqrt{|x|}} \, dx + \int_{-1}^1 \frac{1}{\sqrt{|\sin x|}} - \frac{1}{\sqrt{|x|}} \, dx$$

The first integral I can do symbolically, and the second integrand is now well-behaved at x=0, and you should have no difficulties numerically integrating it. Of course, you might have some numerical stability issues computing values of the integrand near x=0, so you'll have to do some more work to figure out how to evaluate it with little error.

 

[D H]

Another way to address this problem is to represent the gravitation potential via spherical harmonics. http://www.cdeagle.com/pdf/gravity.pdf" (200x200, extensible to 360x360). JPL has developed a 150x150 spherical harmonic model of the Moon's gravity field.

 

[Ulysees]

Thanks.

Can't open that article, D H, but can spherical harmonics be generalized for any shape defined in terms of cubic surface patches? This is what I am after eventually.

 

[Ulysees]

$$\int_{-1}^1 \frac{1}{\sqrt{|\sin x|}} \, dx = \int_{-1}^1 \frac{1}{\sqrt{|x|}} \, dx + \int_{-1}^1 \frac{1}{\sqrt{|\sin x|}} - \frac{1}{\sqrt{|x|}} \, dx$$

The first integral I can do symbolically, and the second integrand is now well-behaved at x=0, and you should have no difficulties numerically integrating it.

How can I tell the integral is well-behaved?

 

[Ulysees]

I'm assuming you are treating your original mass as a continuous matter distribution?

Yes.

You can do some actual analysis and estimate errors. Suppose you are decomposing your mass into little cubes. Then, for each cube, you can compute the smallest and largest densities in that cube. Then, for each component of gravitational force, you should be able to determine the largest and smallest forces that are possible for the cube to exert on your test particle under all possible matter distributions satisfying that density constraint. This way, you have actual hard upper and lower bounds for the value you're calculating.

I've done something in that line when I wanted to decide how small to make the integration step given the distance to the point of interest. Ie if the cube is too close, the point mass approximation is no good - increase resolution and repeat with smaller cubes.

 

[D H]

Thanks.

Can't open that article ...

I edited the post. Thanks for reporting the problem.