If the gravitational force were inversely proportional to distance (rather than distance squared), will celestial bodies fall into each other? If gravity was inversely proportional to distance, will the dynamics of celestial bodies be much different from our world? Will celestial bodies fall into each other?
 A: I think there are two ways to interpret your question - joseph h has given an excellent explanation of what were to happen if you had the force of gravity proportional to $1/r$ instead of the physical $1/r^2$.
My first thought when I read your question is that you mean to ask "what if gravity had the opposite relation to distance than what we observe?" In that case, we would have something like
$$
F_g = GMmr^2
$$
where now $G$ has different units than usual. This means that objects further apart would experience a stronger pull due to gravity than objects close together.
The speed of an object orbiting another would be given by
$$
v = \sqrt{GMr^3}.
$$
This would be madness! It states that the larger an orbit, the faster the object must go in order for that orbit to be stable. Objects very close to the mass they're orbiting would move slower. Neptune would have an orbital speed of $10^{29}$ ms$^{-1}$, far in excess of the speed of light.
Let's be glad that gravitational force decreases with distance, and not the inverse.
Now, would massive objects fall into each other? I think it's hard to say without running a simulation. My intuition leads me to this, though:
Consider a universe with only two objects, placed far apart. At first the attraction is very strong and they move towards each other. Over time, the force will become weaker, and assuming there is nothing to change their velocity, so they would collide. If we give one of the objects an initial sideways velocity, they would not collide, because as they move closer together, the force pulling them together is weaker, so they would likely initially miss each other, and then be attracted again as their distance increased. I think, over a very large time scale, the objects would settle into an equilibrium state where they are at rest with zero distance between them.
A: Why not test it? The following Mathematica code numerically integrates the equations of motion for $F\propto 1/r$.
G = 1; M = 1;
T = 20;
r0 = 1;
dv = .1;
sols = NDSolve[{ x''[t] == -((G M)/(x[t]^2 + y[t]^2)) x[t], 
    y''[t] == -((G M)/(x[t]^2 + y[t]^2)) y[t], x[0] == r0, y[0] == 0, 
    x'[0] == 0, y'[0] == Sqrt[G M] + dv}, {x, y}, {t, 0, T}];
ParametricPlot[
 Evaluate[{{Cos[t], Sin[t]}, {x[t], y[t]}} /. sols], {t, 0, T}, 
 AspectRatio -> 1]

where $T$ is the integration time, $r0$ is the starting radius and $dv$ is the deviation from a circular trajectory. The circular trajectory has been calculated with the help of @joseph h’s answer. This code gives the following plots for different $dv$:

The blue circle shows the reference circular trajectory. We notice two important things. Firstly the orbits precess. They generally don't end up at their starting point. Non-precessing orbits are a special characteristic of Keplerian orbits. Secondly the orbits are still bound. They don't spiral inward. To make sense of this we can look at the potential $V(r)=GM\log(r)$. If you plot this and compare it to $V(r)=-\frac{GM}{r}$ they actually have very similar shape. But, because $\log(r)$ doesn't have an asymptote, there are no longer escape trajectories. Every orbit will eventually return even though it will take very long to return for large velocities.
Edit: for reference I will also include these plots for a $F\propto 1/r^2$ to get a sense of how large $dv$ is. At $dv=+0.5$ we already have an escape trajectory.

A: This is not a complete answer but I feel it is a contribution. The n-body problem is a standard candidate for parallel simulations on GPUs, so it is relatively easy to look into.
Now Apple releases this source to show off its Metal framework. In the file AAPLKernels.metal we find the following function:
static float3 computeAcceleration(const float4 vsPosition,
                                  const float4 oldPosition,
                                  const float  softeningSqr)
{
    float3 r = vsPosition.xyz - oldPosition.xyz;

    float distSqr = distance_squared(vsPosition.xyz, oldPosition.xyz);

    distSqr += softeningSqr;

    float invDist  = rsqrt(distSqr);
    float invDist3 = invDist * invDist * invDist;

    float s = vsPosition.w * invDist3;

    return r * s;
}

It basically calculates the acceleration as $\vec{r} \frac{m}{|r|^3}$ (with the addition of softeningSqr for numerical reasons). I have modified the function to calculate instead $\vec{r} \frac{m}{|r|^2}$ just modifying two rows:
float invDist2 = invDist * invDist;

float s = vsPosition.w * invDist2;

The result of the simulation when the force is inversely proportional to the distance squared:

When I modify the function to use a force inversely proportional to the distance:

Now the $1/r$ shows a different dynamics, like the particle 'bounce off' each other more than in the $1/r^2$ case and coagulate differently. I don't know the reason, but I suggest you, if you have a good GPU, to take a look at this website and modify the code of the shader in the same way. Unfortunately the task is too heavy for my 6-year-old no-GPU laptop.
A: The best way to think about weird modifications to gravity is via the effective potential.  Basically, it is always possible to re-cast the radial motion of a particle experiencing a central force as a 1-D problem of a particle moving in a potential
$$
U_\text{eff}(r) = \frac{l^2}{2m r^2} + U(r),
$$
where $l$ is the angular momentum of the particle, $m$ is its mass, and $U(r)$ is the "real" potential energy of the particle.
So suppose that the force between the particle and the center was $F = k/r$ for some $k$.  This implies a potential of $U(r) = k \ln(r/r_0)$ for some reference radius $r_0$.  (Changing the value of $r_0$ just changes $U(r)$ by a constant, which you'll recall does not make a physical difference.)  This means that the effective potential is
$$
U_\text{eff}(r) = k \ln \left(\frac{r}{r_0}\right) + \frac{l^2}{2 m r^2},
$$
which will look something like this for any non-zero values of the constants involved:

The radial motion of this particle (in 3-D) will be exactly analogous to a ball rolling around in this 1-D potential.  In particular, we can see two things:

*

*The particle can never get to $r = 0$ so long as $l \neq 0$.  (The "real" potential goes to $-\infty$ as $r \to 0$, but the $l^2/2mr^2$ piece of the effective potential goes to $+\infty$, and it's not hard to show that the positive piece always wins out.)


*The particle can never get out to $r \to \infty$, since $U(r) \to \infty$ and $r \to \infty$.
This implies that for any value of the particle's energy, it will oscillate between some minimum and maximum radius indefinitely.  So long as it has a non-zero angular momentum, it cannot spiral in to the center, and it cannot escape to infinity.
The particle will also be going around the center (in the tangential direction) at the same time.  It can be shown, however, that most such orbits will not form closed curves.  Instead, most of the orbital curves will be space-filling, eventually getting arbitrarily close to any point in the orbital plane that lies between the minimum and maximum allowable distances.
A: The classical description of gravitational force is given by Newton’s law of universal gravitation
$$F =\frac{GMm}{r^2}$$
where $F$ is the force between two bodies of mass $M$ and $m$ and $r$ is the distance between their centres and $G$ is the gravitational constant.
So gravitational force is proportional to $1/r^2$. This would mean that if you doubled the distance of one object, the force would reduce by a factor of a quarter. Now if it were the case that*
$$F \stackrel{?}{=} \frac{GMm}{r}$$
(and analogous changes were made to general relativity) then doubling the distance of one object would result in the force decreasing by a half. This seems to suggest that orbiting objects would spiral closer since the force is now increased by a factor of $r$ from what it was. And let’s consider what would happen to the speed of objects in orbit*:
$$\frac{GMm}{r} \stackrel{?}{=} \frac{mv^2}{r}$$
meaning
$$v \stackrel{?}{=} \sqrt{GM}$$
when it used to be inversely proportional to $\sqrt{r}$ in the Newtonian case. This would mean that there would be no dependence on the distance the orbiting object is from the larger mass $M$. This would seem to suggest that orbits would be unstable.
Considering this and what was discussed above, how the gravitational force is now increased by a factor of $r$, you may be correct when you say celestial bodies will fall toward each other, and perhaps we would see this effect of “gravitational bodies falling into each other”, as you put it.
At the very least the universe would be substantially different to what it is now, assuming this new law would allow the universe to exist for any substantial amount of time to begin with.

*

*These equations obviously have incorrect dimensions too.

