Why isn't the work done by gravity positive in this situation? I want to find the work done by the force of gravity to move an object of mass $m$ from infinity to a point $P$ at distance $r_p$ from a body of mass $M$ (that I assume fixed). The formula should be
\begin{equation*}
W = \int_{\infty}^{P} \vec{F}(r) \cdot d\vec{r}
\end{equation*}
The force is conservative, so I can assume a straight path, and since the force on $m$ is directed toward the mass $M$, it is parallel to the displacement and the dot product is positive, equal to magnitude of the force.
\begin{equation*}
\vec{F}(r) \cdot d\vec{r} = G\frac{Mm}{r^2}dr
\end{equation*}
But doing so result in a negative work
\begin{equation*}
\int_{\infty}^{P} G\frac{Mm}{r^2} dr = -G\frac{mM}{r_p}
\end{equation*}
and this is wrong since it should be positive. What accounts for the discrepancy?
Note that Halliday and Resnick in their book "Physics" (chapter 13), do almost the same thing, because they compute the work from the point $P$ to infinity. In their case the math checks out because the force is antiparallel to the displacement and the dot product is negative, resulting in a negative work (which is ok because in that case the movement is against the force). Here it should be positive, but the problem is that not only the integral bounds are switched, but also the sign of the integrand, and so the overall sign is not changed.
 A: This is kind of a weird one, so let's follow the logic step by step:
First of all, you're right that $\vec{F}\cdot\mathrm{d}\vec{r}$ should be positive. It has to be, because as you noticed, the force is always acting in the same direction as the path of integration, so each infinitesimal contribution to the integral is going to be positive. If you expressed it as a Riemann sum, you'd be summing up a bunch of positive quantities.
Given that $\vec{F}\cdot\mathrm{d}\vec{r} > 0$, and $\vec{F}$ obviously points in the negative $\hat{r}$ direction, it must also be the case that $\mathrm{d}\vec{r}$ points in the negative $\hat{r}$ direction. You might say
$$\mathrm{d}\vec{r} = -\lvert\mathrm{d}r\rvert\hat{r}$$
But here's the tricky part. That last equation actually involves two different variables. The $\mathrm{d}\vec{r}$ on the left side is a differential that represents an infinitesimal progression along the path of integration (what might be otherwise denoted $\mathrm{d}\vec{\ell}$ or $\mathrm{d}\vec{s}$), while the $r$ on the right side is a coordinate which measures distance from the origin (what might otherwise be denoted $\rho$, or $x$ if you're integrating along the $x$ axis, or so on).

So to avoid confusion, you might want to write
$$\mathrm{d}\vec{\ell} = -\lvert\mathrm{d}\rho\rvert\hat{\rho}$$
Now, as you step along the line of integration from $\infty$ towards $P$, $\rho$ decreases. So if you're going to do a change of variables to replace $\mathrm{d}\vec{\ell}$ with its expression in terms of $\rho$, you're going to need $\mathrm{d}\rho < 0$, which in turn means that $\mathrm{d}\vec{\ell} = \mathrm{d}\rho\,\hat{\rho}$. Therefore,
$$\int_{\infty}^{P}\vec{F}\cdot\mathrm{d}\vec{\ell} = \int_{\rho=\infty}^{\rho=r_P}\biggl(-G\frac{Mm}{\rho^2}\hat{\rho}\biggr)\cdot(\mathrm{d}\rho\,\hat{\rho})$$
The difference with what you were trying to do was that you assumed $\mathrm{d}r$ (non-vector, after the dot product) would be positive. That $\mathrm{d}r$ actually corresponds to $\mathrm{d}\rho$, and is negative. Roughly speaking, you implicitly did a change of variables from $\ell$ to $\rho$, which introduces a sign change, but it was hidden from you because you used the same letter ($r$) for both variables.
A precise formulation
Here's a fully precise way to do this: take the definition of work as a line integral of the force vector field,
$$W = \int_{0}^{1} \vec{F}\bigl(\vec{r}(s)\bigr)\cdot\underbrace{\vec{r}'(s)\,\mathrm{d}s}_{\mathrm{d}\vec{\ell}}$$
In this definition, $s$ is a parameter of the path of integration, which ranges from $0$ at the beginning of the path to $1$ at the end. So in your case, $s = 0$ at infinity and $s = 1$ at $P$. One possible realization of this, if you assume that $P$ is a point on the $x$ axis, is
$$\vec{r}(s) = \biggl(\frac{r_P}{s}, 0, 0\biggr)$$
Using the above realization, we have
$$\begin{align}
\vec{r}'(s) &= \biggl(\underbrace{-\frac{r_P}{s^2}}_{\mathrm{d}\rho/\mathrm{d}s}, 0, 0\biggr) = -\frac{r_P}{s^2}\hat{x} &
\vec{F}\bigl(\vec{r}(s)\bigr) &= -G\frac{Mm}{(r_P/s)^2}\hat{x}
\end{align}$$
Note that this makes the direction of $\vec{F}$ radially inward, as it should be, at all points along the line of integration.
Plugging these into the formula for work, we get
$$W = \int_{0}^{1} \biggl(-G\frac{Mm}{(r_P/s)^2}\hat{x}\biggr)\cdot\biggl(-\frac{r_P}{s^2}\hat{x}\biggr)\,\mathrm{d}s = G\frac{Mm}{r_P}\int_0^1\mathrm{d}s = G\frac{Mm}{r_P}$$
"Guide questions" for self-study
In case you'd like to ponder this more deeply, here are some "guide questions" to think about:

*

*Convince yourself that, as $s$ ranges from $0$ to $1$, $\vec{r}(s)$ traces out the path on which you want to integrate.

*How does each component of $\vec{r}$ vary as $s$ changes from $0$ to $1$? Larger, smaller, or the same? Try making a drawing that shows it.

*What are some other possible formulas for $\vec{r}(s)$ that would describe the same path?

*What are some other paths, and the formulas that describe them, which would also be suitable for this physical situation? (E.g. if you don't assume $P$ is on the $x$ axis)

*What is the sign of the x-component of $\vec{r}'(s)$?

*Is that consistent with how your position changes as you move along the integration path?

A comparison
It's been pointed out that an answer to another, similar question says that $\mathrm{d}\vec{r}$ is not a vector directed along the path. After looking into it, I believe that's not correct, for a similar reason to what you did here: that answer uses the definition
$$\mathrm{d}\vec{r} = \mathrm{d}r\,\hat{r} + r\,\mathrm{d}\theta\,\hat{\theta} + r\sin\theta\,\mathrm{d}\phi\,\hat{\phi}$$
but that expression is for the differential of the coordinate vector (of which my $\mathrm{d}\rho\,\hat{\rho}$ is the radial part). $\mathrm{d}\vec{\ell}$ is a different entity, which is directed along the integration path.
A: $\begin{equation*}
W = \int_{\infty}^{P} \vec{F}(r) \cdot d\vec{r}
\end{equation*}$
$\vec r = r \,\hat r$ is the displacement from the centre of the mass $M$ and $d\vec r=dr\,\hat r$ is incremental displacement.
The gravitational force is attractive and so in the opposite direction to $\hat r$ thus $\vec F(r) = -G \dfrac {Mm}{r^2} \hat r$
The integral now becomes $\displaystyle \int _\infty^r \left (-G \dfrac {Mm}{r^2} \hat r\right)\cdot (dr\,\hat r)\int _\infty^r -G \dfrac {Mm}{r^2}\,dr$ which will give you your required positive value.
