Singularity Handling in Numerical Integration I'm working on a school project where I'm trying to find stable orbits around a homogenous toroidal planet. The toroid in question has major radius R=2 and minor radius r=1. The density, gravitational constant, and mass of the orbiting object are not really important, I'm just having trouble with the integration.
I've represented the toroid with these parametric equations:
$$x=(2+\rho*cos(\phi))cos(\theta)$$
$$y=(2+\rho*cos(\phi))sin(\theta)$$
$$z=\rho*sin(\phi)$$
where:
$$0\le\theta\le2\pi$$
$$0\le\phi\le2\pi$$
$$0\le\rho\le1$$
I found the jacobian without any problems, and I can integrate over the region fine. I want to find the gravity vector from the toroid at any given test point, and I'm doing each component individually, starting with x, then y, then z. This works fine for far away points, but within the "hole of the donut" there is a singularity whenever the distance from the test point to the area being integrated is zero.
For example: finding the x-component of the gravitational acceleration vector at the test point {0,0,0}
$$\int_0^{2\pi}\!\!\!\int_0^{2\pi}\!\!\!\int_0^1\frac{1}{|x|*x}dV$$
should evaluate to zero. Mathematica can't evaluate it absolutely, and numerically it chokes on the singularity where x=0. I know this integral diverges:
$$\int_0^1\frac{1}{x^2}dx$$
but this integral is zero:
$$\int_{-1}^1\frac{1}{x^2}dx$$
so I'm pretty sure I should get an answer for my gravity integral. What I'm asking is how can I handle the singularity?
If it helps, here is the mathematica input and output:
Ssqr[x_] := Abs[x]*x
x[u_, v_, w_] := (2 + u*Cos[w])*Cos[v]
NIntegrate[1/Ssqr[x[u, v, w]]*jac, {u, 0, 1}, {v, 0, 2 \[Pi]}, {w, 0, 2 \[Pi]}]

(* 
  Out: NIntegrate::slwcon: *slow convergence message*
  2.66615*10^16
*)

I also tried the following with no luck:
NIntegrate[1/Ssqr[x[u, v, w]]*jac, 
  {u, 0, 1}, {v, 0, 2 Pi}, {w, 0, 2 Pi}, 
  Exclusions -> x[u, v, w] == 0]

 A: If I understand your problem correctly, I'd set it up something like so:
x[rho_, phi_, theta_] = (2 + rho*Cos[phi]) Cos[theta];
y[rho_, phi_, theta_] = (2 + rho*Cos[phi]) Sin[theta];
z[rho_, phi_, theta_] = rho*Sin[phi];
p[rho_, phi_, theta_] = {
   x[rho, phi, theta],
   y[rho, phi, theta],
   z[rho, phi, theta]
   };
j[rho_, phi_, theta_] = Simplify[Det[
    {Derivative[1, 0, 0][p][rho, phi, theta],
     Derivative[0, 0, 1][p][rho, phi, theta],
     Derivative[0, 1, 0][p][rho, phi, theta]}
    ]];
j[\[Rho], \[CurlyPhi], \[Theta]] // TeXForm

$$\rho  (\rho  \cos (\varphi )+2)$$
Thus, your integrand should be
j[rho, phi, theta]/(Abs[x[rho, phi, theta]] x[rho, phi, theta]) // TeXForm

$$\frac{\rho  \sec (\theta )}{\left| (\rho  \cos
   (\varphi )+2) \cos (\theta )\right| }$$
So we're looking at
$$\int_0^{2\pi}\int_0^{2\pi}\int_0^1 
\frac{\rho  \sec (\theta )}{\left| (\rho  \cos
   (\varphi )+2) \cos (\theta )\right| } d\rho \, d\varphi \, d\theta$$
which can be separated into
$$\left(\int_0^{2\pi}\int_0^1 
\frac{\rho}{\left| (\rho  \cos
   (\varphi )+2)\right| } d\rho \, d\varphi \right)
\left(\int_0^{2\pi} 
\frac{\sec(\theta )}{\left|\cos(\theta )\right| } d\theta \right).$$
Now, the first integral in the product can be handled analytically
Integrate[rho/Abs[rho*Cos[phi] + 2], {rho, 0, 1}, {phi, 0, 2 Pi}] // TeXForm

$$-2 \left(\sqrt{3}-2\right) \pi$$
The second one cannot be done analytically, at least not by Mathematica's Integrate command, but can be done numerically.
NIntegrate[
  Sec[theta]/Abs[Cos[theta]],
  {theta, 0, 2 Pi},
  Method -> "PrincipalValue",
  Exclusions -> {Pi/2, 3 Pi/2},
  AccuracyGoal -> MachinePrecision
]

(* Out: -4.44089*10^-16 *)

Note that we must specify Method -> "PrincipalValue" together with associated Exclusions.  The specification of the AccuracyGoal effectively overrides the use of PrecisionGoal and is necessary since zero effectively has infinite precision.
