（離散的な）質点間にはたらく万有引力の法則が，連続的な質量密度分布の場合にはどうなるか，特に，質量密度分布が球対称な場合はどうなるか，という話。
ここでは4元ベクトルの出番はなく，すべて3次元ベクトル。なので，太字でベクトルを表し，
$$\boldsymbol{F} = (F_x, F_y, F_z)$$などと書くことにする。
ポアソン方程式までいかずに導く例
原点においた質点がつくる万有引力・重力加速度ベクトル
原点 \(\boldsymbol{r} = \boldsymbol{0}\) においた質量 \(M\) の質点がつくる万有引力 \(\boldsymbol{F}\) を，位置ベクトル \(\boldsymbol{r}\) の位置にある質量 \(m\) の質点が受けるとすると，\( \boldsymbol{F}\) は
$$ \boldsymbol{F} = -\frac{GMm}{r^3} \boldsymbol{r} \equiv m \boldsymbol{g}$$つまり，原点においた質量 \(M\) の質点がつくる重力加速度ベクトル \(\boldsymbol{g}\) は
$$ \boldsymbol{g} = -\frac{GM}{r^3} \boldsymbol{r}$$
任意の位置においた質点がつくる重力加速度ベクトル
位置 \(\boldsymbol{r}_1\) においた質量 \(M_1\) の質点が位置 \(\boldsymbol{r}\) につくる重力加速度ベクトルは，質点までの位置が \(\boldsymbol{r} \rightarrow \boldsymbol{r} -\boldsymbol{r}_1\) に変わるので，
$$ \boldsymbol{g} = -\frac{GM_1}{|\boldsymbol{r} -\boldsymbol{r}_1|^3} \,(\boldsymbol{r} -\boldsymbol{r}_1)$$
2個の質点がつくる重力加速度ベクトル
位置 \(\boldsymbol{r}_1\) においた質量 \(M_1\) の質点と，位置 \(\boldsymbol{r}_2\) においた質量 \(M_2\) の質点が位置 \(\boldsymbol{r}\) につくる重力加速度ベクトルは，重ね合わせの原理から
$$ \boldsymbol{g} = -\frac{GM_1}{|\boldsymbol{r} -\boldsymbol{r}_1|^3} \,(\boldsymbol{r} -\boldsymbol{r}_1) -\frac{GM_2}{|\boldsymbol{r} -\boldsymbol{r}_2|^3} \,(\boldsymbol{r} -\boldsymbol{r}_2)$$
複数個の質点がつくる重力加速度ベクトル
位置 \(\boldsymbol{r}_i\) においた質量 \(M_i\) の複数の質点（\(i = 1, 2, 3, \dots\)）が位置 \(\boldsymbol{r}\) につくる重力加速度ベクトルは，
$$ \boldsymbol{g} = -\sum_{i}\frac{GM_i }{|\boldsymbol{r} -\boldsymbol{r}_i|^3} (\boldsymbol{r} -\boldsymbol{r}_i)$$
連続的な質量分布の場合
連続的な質量分布 \(\rho(\boldsymbol{r})\) の場合は，以下のような置き換えをして
\begin{eqnarray}
\boldsymbol{r}_i &\rightarrow& \boldsymbol{r}^{\prime} \\
M_i &\rightarrow& \rho(\boldsymbol{r}_i)\,dV_i \ \ \rightarrow\ \ \rho(\boldsymbol{r}^{\prime})\,dV_i \\
\sum_{i} dV_i &\rightarrow& \iiint dV^{\prime}
\end{eqnarray}
\begin{eqnarray}
\boldsymbol{g} &=& -\sum_{i}\frac{GM_i}{|\boldsymbol{r} -\boldsymbol{r}_i|^3} (\boldsymbol{r} -\boldsymbol{r}_i)\\
&\Rightarrow& -\iiint dV^{\prime} \frac{G \rho(\boldsymbol{r}^{\prime}) }{|\boldsymbol{r} -\boldsymbol{r}^{\prime}|^3}(\boldsymbol{r} -\boldsymbol{r}^{\prime})
\end{eqnarray}
球対称な質量分布がつくる万有引力
原点を中心とした球対称な電荷分布を表す電荷密度は
$$\rho(\boldsymbol{r}) = \rho(\sqrt{x^2 + y^2 + z^2}) \Rightarrow \rho(r)$$
と書ける。
簡単のために，
$$\boldsymbol{r} = (x, y, z) = (0, 0, r)$$とし，\(z\) 軸上の重力加速度ベクトルを求める。一旦求めてしまえば，（球対称性から）\(z\) 軸上に限らず任意の場所での重力加速度ベクトルの答えになる。
また，積分変数を極座標で書くと
\begin{eqnarray}
x’ &=& r’ \sin\theta’ \cos \phi’ \\
y’ &=& r’ \sin\theta’ \sin \phi’\\
z’ &=& r’ \cos\theta’ \\
dV’ &=& (r’)^2 \,dr’\,\sin\theta’ d\theta’ \,d\phi’
\end{eqnarray}
また，
$$|\boldsymbol{r} -\boldsymbol{r}’| = \left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{1}{2}}$$
したがって，
$$\boldsymbol{g} = -\iiint \frac{G \rho(r’) (\boldsymbol{r} -\boldsymbol{r}’)}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}}\, dV’$$
各成分ごとに計算すると，\(x\) 成分は
\begin{eqnarray}
g_x &=& -\iiint \frac{G \rho(r’) (x -x’)}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \,dV’ \\
&=& -\int_0^{\infty} (r’)^2 dr’ \int_0^{\pi} \sin\theta’ d\theta’ {\color{blue}{\int_0^{2\pi} d\phi’}} \frac{G \rho(r’) (0 -r’ \sin\theta’ {\color{blue}{\cos\phi’}})}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \\
&=& 0
\end{eqnarray}
\(y\) 成分も同様にして
\begin{eqnarray}
g_y &=& -\iiint \frac{G \rho(r’) (y -y’)}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \,dV’ \\
&=& -\int_0^{\infty} (r’)^2 dr’ \int_0^{\pi} \sin\theta’ d\theta’ {\color{blue}{\int_0^{2\pi} d\phi’}} \frac{G \rho(r’) (0 -r’ \sin\theta’ {\color{blue}{\sin\phi’}})}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \\
&=& 0
\end{eqnarray}
\(z\) 成分は
\begin{eqnarray}
g_z &=& -\iiint \frac{G \rho(r’) (z -z’)}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \,dV’ \\
&=& -\int_0^{\infty} (r’)^2 dr’ \int_0^{\pi} \sin\theta’ d\theta’ \frac{G\rho(r’) (r-r’ \cos\theta’ )}{\left\{ r^2 + (r’)^2 -2 r r’ \cos\theta’ \right\}^{\frac{3}{2}}} \int_0^{2\pi} d\phi’ \\
&=& -\int_0^{\infty} (r’)^2 dr’ \frac{2\pi G\rho(r’)}{r^2}
\left( \frac{r + r’}{|r + r’|} + \frac{r -r’}{|r -r’|} \right) \\
&=& -\frac{G}{r^2} \int_0^{r} 4\pi \rho(r’) (r’)^2 dr’ \\
&=& -\frac{GM_{\color{blue}{r}}}{r^2} = -\frac{GM_{\color{blue}{r}}}{r^3} r
\end{eqnarray}
ここで \(\displaystyle M_{\color{blue}{r}} \equiv \int_0^{{\color{blue}{r}}} 4\pi \rho(r’) (r’)^2 \,dr’ \) は半径 \({\color{blue}{r}}\) の球内の全質量である。
以上の結果をまとめると，\(\boldsymbol{r} = (0, 0, r)\) として
$$\boldsymbol{g} = -\frac{G M_{\color{blue}{r}} }{r^3}\boldsymbol{r}$$
いったんこのようなベクトル式で表された解は，一般に \(\boldsymbol{r} = (x, y, z)\) としても成り立つことに注意。最終的に，万有引力 \(\boldsymbol{F}\) は
$$\boldsymbol{F} = m \boldsymbol{g} = -\frac{G M_{\color{blue}{r}} m}{r^3}\boldsymbol{r}$$
球対称に分布していれば，半径 \({\color{blue}{r}}\) の外側にどのように分布していても，万有引力 \(\boldsymbol{F}\) は半径 \({\color{blue}{r}}\) 内の全質量 \(M_{\color{blue}{r}} \) で決まり，外部に分布している質量の影響は受けない！というところがおもしろい。ガウスの定理を使わなくても，全空間での積分をきちんとすれば同様の答えが出てくるのです。
この結果は，静電磁場に対するガウスの法則に対応している。以下の電磁気学のページを参照。
静電場との類推から導く例
（離散的な）点電荷に対するクーロンの法則
まずは事前準備として，万有引力ではなく，2つの点電荷（大きさが無視できる質点が電荷を帯びたもの）にはたらく電場による力を考えることからはじめる。
原点においた点電荷 \(Q\) がつくる電場によって，位置ベクトル \(\boldsymbol{r}\) にある点電荷 \(q\) が受ける力 \(\boldsymbol{F}\) は
$$\boldsymbol{F} = \frac{Q q}{4\pi\varepsilon_0} \frac{\boldsymbol{r}}{r^3}$$
連続的な電荷密度分布に対するクーロンの法則
電荷分布が（離散的な）点電荷ではなく，連続分布する電荷密度分布 \(\rho(\boldsymbol{r})\) であらわされる場合は，電荷 \(q\) が受ける力 \(\boldsymbol{F}\) は以下のように表される。
\begin{eqnarray}
\nabla\cdot\boldsymbol{E} &=& \frac{\rho}{\varepsilon_0} \tag{1}\\
\nabla\times\boldsymbol{E} &=& \boldsymbol{0}, \quad\Rightarrow \boldsymbol{E} \equiv -\nabla\phi \tag{2}\\
\therefore\ \ \nabla^2 \phi &=& -\frac{\rho}{\varepsilon_0} \tag{3}\\
\boldsymbol{F} &=& q \boldsymbol{E} \tag{4}
\end{eqnarray}
\((1)\) 式および \((2)\) 式が静電場に対するマクスウェル方程式であり，静電ポテンシャル \(\phi\) に対する \((3)\) 式がポアソン方程式と呼ばれる。
ポアソン方程式の完全な解
$$\phi(\boldsymbol{r}) = \frac{1}{4\pi \varepsilon_0} \iiint \frac{\rho(\boldsymbol{r}’)}{|\boldsymbol{r} -\boldsymbol{r}’|} \,dV’$$を使うと，一般に \((4)\) 式で表される電場による力は
\begin{eqnarray}
\boldsymbol{F} &=& -q \nabla\phi \\
&=& \frac{q}{4\pi \varepsilon_0} \iiint\frac{\rho(\boldsymbol{r}’) (\boldsymbol{r} -\boldsymbol{r}’) }{|\boldsymbol{r} -\boldsymbol{r}’|^3}dV’
\end{eqnarray}
となるのだが…
球対称分布の場合のガウスの定理による解
電荷密度 \(\rho\) の分布が球対称である場合は，\((1)\) 式を半径 \(r\) の球で体積積分し，ガウスの定理を使うと以下のように解ける。
まず，電荷密度分布 \(\rho\) が球対称，つまり動径座標 \(r\) のみの関数であるとすると，\((3)\) 式のポアソン方程式の両辺を見比べると，静電ポテンシャル \(\phi\) もまた球対称，つまり \(r\) のみの関数となるはずである。
$$\rho = \rho(r), \ \ r = |\boldsymbol{r}| \quad\Rightarrow\quad \phi = \phi(r)$$
すると，\((2)\) 式より，電場 \(\boldsymbol{E}\) は動径方向 \(\boldsymbol{r}\) に平行な成分のみを持つベクトルになる。
$$\boldsymbol{E} = -\nabla\phi(r) = -\frac{d\phi}{dr} \nabla r = -\frac{d\phi}{dr} \frac{\boldsymbol{r}}{r} \equiv E_r(r) \frac{\boldsymbol{r}}{r} $$
次に，\((1)\) 式の両辺を原点を中心とした半径 \(r\) の球の体積 \(V\) で体積積分する。
$$\iiint_V (\nabla\cdot\boldsymbol{E}) \,dV = \frac{1}{\varepsilon_o} \iiint_V \rho(r) \,dV$$
右辺の積分は，以下のようにして半径 \(r\) の球内にある全電荷 \(Q_r\) で書ける。
$$ \iiint_V \rho(r) \,dV = 4\pi \int_0^r \rho(r) r^2 dr \equiv Q_r$$
左辺はガウスの定理により，体積 \(V\) の表面である，半径 \(r\) の球面 \(S\) での面積積分に帰着して…
\begin{eqnarray}
\iiint_V (\nabla\cdot\boldsymbol{E})\, dV &=& \iint_S \boldsymbol{E}\cdot\boldsymbol{n}\, dS
\qquad\left(\boldsymbol{n} = \frac{\boldsymbol{r}}{r}\right) \\
&=& \iint_S E_r(r) \,dS \\
&=& 4 \pi r^2 E_r(r)
\end{eqnarray}
両辺を等しいとおいて
\begin{eqnarray}
4 \pi r^2 E_r(r) &=& \frac{Q_r}{\varepsilon_0} \\
\therefore\ \ E_r(r) &=& \frac{Q_r}{4 \pi \varepsilon_0 r^2 }\\
\therefore\ \ \boldsymbol{E} &=& E_r(r) \frac{\boldsymbol{r}}{r} = \frac{Q_{r}}{4\pi \varepsilon_0}\frac{\boldsymbol{r}}{r^3}
\end{eqnarray}
したがって，この場合のクーロンの法則は
$$\boldsymbol{F}(\boldsymbol{r}) = \frac{Q_r q}{4\pi\varepsilon_0} \frac{\boldsymbol{r}}{r^3}$$
となり，原点に点電荷 \(Q\) をおいた場合のクーロンの法則で \(Q\) を \(Q_r\) に置き換えた形になることに注意。
（離散的な）質点に対する万有引力の法則
2つの質点にはたらく万有引力を考える。
原点においた質量 \(M\) の質点（大きさが無視できる粒子状？の物体）によって，位置ベクトル \(\boldsymbol{r}\) にある質量 \(m\) の質点が受ける力 \(\boldsymbol{F}\) は
$$\boldsymbol{F} = \frac{GMm}{r^3} \boldsymbol{r}$$
これが，物体 $M$ と $m$ が空間を隔てて直接，力を及ぼすという遠隔作用としての万有引力の法則である。物体 $M$ と $m$ の間の空間には何もないし，物体 $M$ が及ぼす力は，どれだけ離れていても物体 $m$ に瞬時に伝わる，とする。
連続的な質量密度分布に対する万有引力の法則
連続的な質量密度分布 \(\rho(\boldsymbol{r})\) であらわされる場合は，上記の電磁気学的扱いをそのままパクって
\begin{eqnarray}
\nabla\cdot\boldsymbol{g} &=& -4 \pi G \rho \\
\boldsymbol{g} &\equiv& -\nabla\phi \\
\therefore\ \ \nabla^2 \phi &=& 4\pi G \rho\\
\boldsymbol{F} &=& m \boldsymbol{g}
\end{eqnarray}
ここで \(\boldsymbol{g}\) は（あえて名前をつけるとすれば）重力加速度ベクトルである。
これらの方程式が，質量密度分布 \(\rho(\boldsymbol{r})\) が周辺に，まず（不均一な）重力ポテンシャル場 $\phi$ をつくり，その場の勾配が力として作用するという，場の理論，近接作用としての万有引力の法則である。
物体 $M$ が空間に存在すると，その周辺一帯には（不均一な）重力ポテンシャル場ができる。その空間のどこかに物体 $m$ が置かれると，置かれた場所の重力ポテンシャル場の不均一さを感知して，勾配にそって運動する。物体 $m$ は離れた場所にある物体 $M$ からの作用を直接受けるのではなく，自身の周辺につくられた重力ポテンシャル場の不均一さを感知するのだという意味で，遠隔作用とは異なる。（ただし，物体 $M$ がつくる重力ポテンシャル場は瞬時にして全空間に広がるとしか解釈できないので，結果的には，瞬時にして物体 $M$ の作用を受けてしまうのであり，電磁場のような本格的？な古典場・近接作用場という解釈までにはいかないところが何とも難しい。）
電磁気学的類推により，質量密度 \(\rho\) の分布が球対称である場合の万有引力 \(\boldsymbol{F}\) は
$$\boldsymbol{F}(\boldsymbol{r}) = \frac{G M_r m}{r^3} \boldsymbol{r}$$
となり，原点に質点 \(M\) をおいた場合の万有引力の法則で， \(M\) を半径 \(r\) の球の体積 \(V\) 内の全質量 \(M_r\)
$$M_r \equiv \iiint_V \rho(r) \,dV = 4\pi \int_0^r \rho(r) r^2 dr$$
に置き換えた形になることは明らかであろう。