Minimizing a potential using Solve

Question

I am having problems minimizing a potential:

$\text{V}(h,\eta)=\gamma \left(-h^2\right) \left(\eta ^2 \cos ^2(\theta )+\eta \cos (\delta ) \sqrt{-\eta ^2-h^2+1} \sin (2 \theta )+\left(-\eta ^2-h^2+1\right) \sin ^2(\theta )\right)$ (Input code below)

V = -h^2 γ (η^2 Cos[θ]^2 + (1 - h^2 - η^2) Sin[θ]^2 + η Sqrt[1 - h^2 - η^2] Cos[δ] Sin[2 θ])

I try to minimize simply using:

sol = Solve[{D[V, h] == 0, D[V, η] == 0}, {h, η}]

It seems to solve quickly and with no problems giving me 9 solutions, these solutions depend on the parameters $\theta, \delta, \gamma$.

However if for example I try to calculate:

Chop[N[D[V,h] /. sol /. θ -> 1 /. δ -> 1 /. γ -> 1]]

The output it gives is:

{0, -0.104069, 0, 0, -0.484451, 0.104069, 0, 0, 0.484451}

As you can see extrema 2, 5, 6 and 9 are not even close to zero! It seems Mathematica is solving this incorrectly, or maybe the solutions are only valid in specific regions (not the full domain)?

Does anyone have any ideas? I thought maybe the problem was with the square root in the function (otherwise it's a simple polynomial) but have tried solving it using the Lagrange multiplier method taking a function $V(h,η,s)$ and $s=\sqrt{1-h^2-η^2}$ but it also doesn't seem to give me all zero solutions when I sub it in to the Lagrange equation $\frac{\partial V}{\partial h}+\lambda \frac{\partial g}{\partial g}$, where $g=s-\sqrt{1-h^2-η^2}=0$, the constraint applied to this function.

Any ideas what's happening ? Has anyone encountered this problem before?

Answer 1

There are minimum example of your problem

f = Sqrt[a^2 - x^2] + x;
Plot[{f /. a -> 1}, {x, -1, 1}]

It obviously has only one extremum, but Solve returns two

sol = Solve[D[f, x] == 0, {x}]
D[f, x] /. sol /. a -> 1

{{x -> -(a/Sqrt[2])}, {x -> a/Sqrt[2]}}
{2, 0}

It is because Sqrt is a multivalued function in the complex domain. Let's plot both Riemann surface sheets of Sqrt

f2 = -Sqrt[a^2 - x^2] + x;
Plot[{f /. a -> 1, f2 /. a -> 1}, {x, -1, 1}]

Now it is clear that the parasitic solution came from the -Sqrt.

If you solve over Reals domain the result is correct

sol = Solve[D[f, x] == 0, {x}, Reals]
D[f, x] /. sol /. a -> 1

{{x -> Sqrt[a^2]/Sqrt[2]}}
{0}

The original problem is too complicated to solve it over the real domain. But you can always select right solutions manually.

Another possibility is setting parameters θ,δ,γ before Solve

θ = 1;
δ = 1;
γ = 1;
sol = Solve[D[V, h] == 0 && D[V, \[Eta]] == 0, {h, \[Eta]}];
Chop[N[D[V, h] /. sol]]

{0, 0, 0, 0, 0}