Compute the Legendre transform for a singular Lagrangian

As a quick note, the equations of motion that come from that Lagrangian are

$$\frac{d}{dt}\left(\dot q_1 + \dot q_2\right) = -2kq_1^3$$ $$\frac{d}{dt}\left(\dot q_2 + \dot q_1\right) = -2kq_2^3$$

The Dirac procedure for singular Lagrangians goes as follows:

Step 1: Calculate the generalized momenta as usual $$p_1 \equiv \frac{\partial L}{\partial \dot q_1} = \dot q_1+\dot q_2$$ $$p_2 \equiv \frac{\partial L}{\partial \dot q_2} = \dot q_1+\dot q_2 = p_1$$

Clearly this is not invertible. We have one "good" equation (defining $p_1$ in terms of the generalized velocities) and one "bad" equation ($G\equiv p_2-p_1 = 0$, an algebraic relation between the momenta in which the velocities do not appear).

$G=0$ is called a primary constraint of the Dirac procedure - an algebraic relation between momenta and (possibly) coordinates, in which the generalized velocities are absent.

Step 2: Compute the naive Hamiltonian

If we compute the Hamiltonian as usual, we find

$$H_0 = p_1\dot q_1 + p_2 \dot q_2 - L = p_1(p_1-\dot q_2) + p_2 \dot q_2 - \frac{1}{2}p_1^2 + \frac{k}{2}(q_1^4+q_2^4)$$ $$ = \frac{p_1^2}{2} + (p_2-p_1)\dot q_2 + \frac{k}{2}(q_1^4+q_2^4)$$ $$ = \frac{p_1^2}{2} + \frac{k}{2}(q_1^4+q_2^4)$$

If you calculate the Hamilton equations, you will see that they don't match with the Lagrangian equations:

$$\dot p_1 = -\frac{\partial H_0}{\partial q_1} = -2kq_1^3$$ $$\dot p_2 = -\frac{\partial H_0}{\partial q_2} = -2kq_2^3$$ $$\dot q_1 = \frac{\partial H_0}{\partial p_1} = p_1$$ $$\dot q_2 = \frac{\partial H_0}{\partial p_2} = 0$$

Step 3: Extend the phase space and construct the complete Hamiltonian

We now extend the phase space by introducing a new variable $v$, and defining it to Poisson-commute with the regular phase space variables, i.e. $$\{v,q_i\} = \{v,p_i\} = 0$$

The complete Hamiltonian is obtained by multiplying $v$ by our primary constraint $G$ and adding it to $H_0$:

$$H = \frac{p_1^2}{2} + \frac{k}{2}(q_1^4+q_2^4) + v(p_2-p_1)$$ The new Hamiltonian equations are

$$\dot p_1 = -\frac{\partial H}{\partial q_1} = -2kq_1^3$$ $$\dot p_2 = -\frac{\partial H}{\partial q_2} = -2kq_2^3$$ $$\dot q_1 = \frac{\partial H}{\partial p_1} = p_1-v$$ $$\dot q_2 = \frac{\partial H}{\partial p_2} = v$$

Step 4: Obtain additional algebraic relations

Because $G$ is identically zero, it must be that $$\dot G = \dot p_2 - \dot p_1 = 0$$ $$\implies T\equiv q_2^3-q_1^3 = 0$$

We call $T$ a secondary constraint of the Dirac procedure - a constraint obtained through differentiation of a primary constraint, and then simplified by using the Hamilton equations obtained from the complete Hamiltonian (though in this case, the naive Hamiltonian would have done just as well).

Step 5: Determine $v$ and eliminate it from the complete Hamiltonian

Differentiating the secondary constraint allows us to determine $v$:

$$\dot T = 3(q_1^2 \dot q_1 - q_2^2 \dot q_2) = 3(q_1^2[p_1-v] - q_2^2[v])$$ $$= 3(q_1^2 p_1 - (q_1^2+q_2^2)v) = 0$$ $$\implies v = \frac{q_1^2 p_1}{q_1^2 + q_2^2}$$

and so

$$H = \frac{p_1^2}{2} + \frac{k}{2}(q_1^4+q_2^4) + \frac{q_1^2 p_1}{q_1^2+q_2^2}(p_2-p_1)$$

And that's it, we're done.

---

You can confirm that this, along with the primary and secondary constraints, reproduces the proper equations of motion:

$$\dot p_1 = -\frac{\partial H}{\partial q_1} = -2kq_1^3 - \frac{2q_1(1-q_1^2)(p_1p_2-p_1^2)}{q_1^2+q_2^2}$$ $$\dot p_2 = -2kq_2^3 +\frac{2q_2q_1^2(p_1p_2-p_1^2)}{q_1^2+q_2^2}$$ $$\dot q_1 = p_1 + \frac{q_1^2}{q_1^2+q_2^2}(p_2-2p_1)$$ $$\dot q_2 = \frac{q_1^2}{q_1^2+q_2^2} p_1$$ $$G \equiv p_2-p_1 = 0$$ $$T\equiv x_2^3-x_1^3 = 0$$

which simplifies to

$$\dot p_1 = \frac{d}{dt}(\dot q_1 + \dot q_2) = -2kq_1^3 = -2kq_2^3$$

---

In summary, singular Lagrangian systems have several common features

1. Defining equations for the generalized momenta yield (some) algebraic equations between phase space variables which do not include generalized velocities, and the system is therefore non-invertible. These equations are called primary constraints, and their derivatives yield secondary constraints
2. The procedure for obtaining the complete Hamiltonian extends the phase space and uses the new variables a bit like Lagrange multipliers to add the primary constraints to the naive Hamiltonian
3. At least some of the "Lagrange multipliers" can be eliminated from the new Hamilton equations by using the primary and secondary constraints, and the resulting system of equations (Hamilton equations + constraints) reproduces the original dynamics
4. This was not included in this example, but any multipliers which remain undetermined at the end of this procedure enter the solutions as arbitrary functions, which wouldn't have been determined by the Lagrangian equations of motion either.

User J. Murray has already given a nice answer. Let us here summarize how the Dirac-Bergmann analysis would proceed in the (possibly conceptionally simpler) coordinates

$$q^{\pm}~:=~q^1\pm q^2, \qquad p_{\pm}~:=~\frac{p_1\pm p_2}{2}.$$

OP's original Lagrangian then reads

$$ L_0~=~ \frac{1}{2} (\dot{q}^+)^2 -V, \qquad V~=~ \frac{k}{16}\left((q^+)^4+(q^-)^4+6(q^+)^2(q^-)^2\right).$$

Primary constraint:

$$p_-~\approx~0.$$

Original Hamiltonian:

$$H_0~=~\frac{1}{2} p_+^2 +V.$$

Consistency check:

$$ 0~\approx~-\dot{p}_-~\approx~\{H_0,p_-\} ~=~\frac{\partial V}{\partial q^-}~=~\frac{k}{4}q^-\left((q^-)^2+3(q^+)^2\right).$$

Secondary constraint:

$$q^-~\approx~0.$$

Result: Hamiltonian: $$H~=~\frac{1}{2} p_+^2 +\frac{k}{16}(q^+)^4~=~\frac{1}{8} (p_1+p_2)^2 +\frac{k}{16}(q^1+q^2)^4 $$ with 2 second-class constraints: $$p_-~\approx~0~\approx~q^-.$$

References:

1. M. Henneaux & C. Teitelboim, Quantization of Gauge Systems, 1994.