Is there a mathematical derivation of potential energy that is *not* rooted in the conservation of energy? For simplicity I'll consider only gravity, but in general this question only applies to conservative forces.
As per my understanding, the way one gets to the equation for gravitational potential energy is this:
$(1)$ $\ddot{r}=-\frac{MG}{r^2}$ is given experimentally
$(2)$ $KE=\frac{1}{2}m\dot{r}^2$ is derived  theoretically and (thought) experimentally 
$(3)$ Energy is conserved
$(4)$ For a set of generic initial conditions, 


*

*Using $(1)$, find $r(t)$ and then derive $\dot{r}(r)$

*Using $(2)$ find $KE(r)$

*Using $(3)$, $PE(r)=E_0-KE(r)$


Is this a correct method?
If the answer is 'no' my question ends here, if not: surely $(3)$ is by definition correct (circularly so), and so the equation for potential energy is arbitrary? I would be interested in seeing a method that does not do this (or is correct at all, if the answer was 'no').
 A: There is no need for any empirical evidence.  This is pure mathematics.  
Step 1:
Assume a force is conservative.  This means that ${\vec \nabla} \times {\vec F} =0$
Step 2:
Then, via Green's theorem, you know that the quantity $\int_{a}^{b}{\vec F}\cdot d{\vec s}$ does not depend on the path you take from a to b.  (equivalently, this integral is zero if the path corresponds to a closed loop)
Step 3:
Then, since the value of that integral is independent of the path, you can then say that the value of $\int_{a}^{b}{\vec F}\cdot d{\vec s}$ depends ONLY on the points a and b, and therefore, we can conceptualize a field that takes a values $V(a)$ and $V(b)$ and that satisfies $\Delta V_{a\rightarrow b} = -\int_{a}^{b}{\vec F}\cdot d{\vec s}$ 
Step 4: 
Since (we are assuming $\vec F$ is the only force in the universe here), ${\vec F} = m{\vec a}$, we have (forgive my abuse of differentials moving the dt over, it's faster than the more rigorous result using the chain rule):
$$\begin{align}
\int {\vec F}\cdot d{\vec s} &= m\int {\vec a}\cdot d{\vec s}\\
&= m\int \frac{d{\vec v}}{dt}\cdot d{\vec s}\\
&= m\int d{\vec v} \cdot \frac{d{\vec s}}{dt}\\
&=m\int {\vec v}\cdot d{\vec v}\\
&=\frac{1}{2}mv_{f}^{2} - \frac{1}{2}mv_{i}^{2}
\end{align}$$
Step 5:
Thus, putting steps 3 and 4 together, we find that
$$\Delta V = - \Delta KE$$
or, as is more commonly written
$$\Delta KE + \Delta PE = 0$$
So, there are no real assumptions or empirical observations necessary.  All it is is calculus and starting with a force that is conservative (the gravitational force you quote DOES satisfy ${\vec \nabla} \times {\vec F} = 0$, which you can check.).  Note that this method lets you derive the potential energy for ANY conservative force without appealing to conservation of energy--you actually PROVE the latter without assuming it!
A: Do you have a problem with force being the derivative of potential energy (up to a minus sign)? Because then you experimentally verify (gravity for example) that $F=-GmM/r^2$ and now $U=\int\frac{dU}{dr}dr=-\int Fdr = GmM\int \frac{dr}{r^2}=-GmM/r$ (up to a constant).
Now of course conservative forces and energy are linked to this, but there is no circular reasoning. (I may have missed your concern though.)
A: The definition of potential energy of a conservative force is not rooted in energy conservation. For instance, the mechanical energy of a system could be shrinking due to e.g. friction but the notion of potential energy could still be well-defined.
In general, the difference $V(B) - V(A)$ in potential energy (associated with some conservative force ${\bf F}$) between positions $A$ and $B$, is defined to be minus the work $W=\int_{A}^{B} {\bf F}\cdot d{\bf r} $ done by the force to get from $A$ to $B$ along any path.
A: The correct definition of change in potential energy is the opposite of the work done by forces internal to a system. Sloppy textbooks authors leave out the bit about having to first choose a system. In its most general form, work is the line integral of a force along a path. Therefore, change in potential energy, and potential energy itself, is also an integral.
A: Although several correct answers were given here I will try to give an alternative answer that uses your equations.
The equation of motion, obtained either by standard Hamiltonian or Lagrangian methods, for a body with kinetic energy given by your expression (2) is
$$m \ddot{r} = - \frac{\partial V}{\partial r}$$
Comparing with your "experimentally given" equation of motion (1) we obtain
$$\frac{\partial V}{\partial r} = \frac{GMm}{r^2}$$
By direct integration
$$V = - \frac{GMm}{r} + V_0$$
The integration constant $V_0$ can be obtained from boundary conditions: the potential energy is zero for infinitely separated bodies $r\rightarrow \infty$. This gives $V_0=0$
