How to find the equillibrium points using Jacobian and Hessian? Given that I have Jacobian and Hessian matrices of three particles interacting with each other in a harmonic trap through Coulomb's law in a 2D plane, how do I find the equilibrium points of them (I know it should be a triangle)? 
 A: With a potential surface $$U=\left(\sum_{k=1}^3A|\mathbf{r}_k|^2\right)+B\left(|\mathbf{r}_1-\mathbf{r}_3|‌​^{-1}+|\mathbf{r}_2-\mathbf{r}_3|^{-1}+|\mathbf{r}_1-\mathbf{r}_2|^{-1}\right),$$ a guess parametrization is to fix the particle coordinates to form the vertices of an equilateral triangle with a radius $R$, ie
$$\mathbf{r}_k=\left(R\text{cos}\left(\frac{2\pi k}{3}\right),R\text{sin}\left(\frac{2\pi k}{3}\right)\right).$$
This reduces the number of dimensions from 6 to 1, at which point you can stick this guess into the potential $U$ and take the derivative with respect to radius $R$ and set equal to zero.
In Mathematica: 
r[k_] := R {Cos[2 \[Pi] k/3], Sin[2 \[Pi] k/3]};
norm[v_] := Sqrt[v.v];
U = Refine[A (norm[r[1]]^2 + norm[r[2]]^2 + norm[r[3]]^2) + B (1/norm[r[1] - r[3]] + 1/norm[r[2] - r[3]] + 1/norm[r[1] - r[2]]), Assumptions -> R > 0];
Solve[D[U, R] == 0, R]

giving $$R= \frac{\sqrt[3]{B}}{\sqrt[3]{2} \sqrt[6]{3} \sqrt[3]{A}}.$$
So the question is, is this correct? It's a critical point of the 1-d potential, which incidentally looks like this:
Plot[U /. {A -> 1, B -> 1}, {R, 0.02, 5}]


But is it also a minimum of the higher-dimensional potential surface? Evaluating the gradient and computing the eigenvalues of the Hessian at the previously-determined minimum yields
{x1, x2, x3} = {{a, b}, {c, d}, {e, f}}
UHighDim = A (norm[x1]^2 + norm[x2]^2 + norm[x3]^2) + B (1/norm[x1 - x3] + 1/norm[x2 - x3] + 1/norm[x1 - x2]);
Refine[(D[UHighDim, {{a, b, c, d, e, f}, 1}] /. Thread[{a, b, c, d, e, f} -> Flatten@{r[1], r[2], r[3]}]) /. R -> B^(1/3)/(2^(1/3) 3^(1/6) A^(1/3)), Assumptions -> {R > 0, A > 0, B > 0}]
Eigenvalues@Refine[(D[UHighDim, {{a, b, c, d, e, f}, 2}] /. Thread[{a, b, c, d, e, f} -> Flatten@{r[1], r[2], r[3]}]) /. R -> B^(1/3)/(2^(1/3) 3^(1/6) A^(1/3)), Assumptions -> {R > 0, A > 0, B > 0}]


{0, 0, 0, 0, 0, 0}
{6 A, 3 A, 3 A, 2 A, 2 A, 0}

Since $\nabla U=\mathbf{0}$ at our previously guessed equilibrium configuration, the guess equilibrium configuration is also a critical point of the higher-dimensional energy surface. In addition, the Hessian is positive semidefinite, which means that the configuration is stable to second order (since everything curves upwards in higher-dimensional space about that point), with the exception of the mode of motion which has eigenvalue zero, which just corresponds to a uniform rotation of the triangle. So our guess is correct.
You can also calculate the harmonic frequencies of the trap from this, if needed.
