Finding the minimum points of a discrete trigonometric equation

Question

Basically I'm trying to solve the equation

$$A\Big(\sin(\theta_a-\theta_{a+1}) + \sin(\theta_a-\theta_{a-1})\Big) + B\Big(n_a^z \sin(\theta_a) - n_a^x \cos(\theta_a)\Big) = 0$$ subjected to the constraint that $$A\Big(\cos(\theta_a-\theta_{a+1}) + \cos(\theta_a-\theta_{a-1})\Big) + B\Big(n_a^x \sin(\theta_a) + n_a^z \cos(\theta_a) \Big) > 0.$$ In addition I would like all of the angles $\theta_a$ to lie on the interval $(-\pi,\pi)$.

Here $A$ and $B$ are some numerical constants and $a$ labels the sites of a lattice. For concreteness let's say that $a=1,2,3,\dots,N$ where $N$ is some finite even number. Furthermore the numbers $n_a^z = 1$ and $n_a^x = 0$ for all sites except $a = N/2$. At the site $a = N/2$ the two numbers can be chosen arbitrarily. For instance $n_{N/2}^z = -1/\sqrt{2}$ and $n_{N/2}^x = -1\sqrt{2}$.

Furthermore I'm using the periodic boundary conditions $\theta_1 = \theta_N = 0$.

To solve this numerically I have tried using the function FindRoot on the first equation. However, this results in that some of the angles I find are not satisfying the second constraint. Is there a nice way of solving the above two equations numerically for a finite lattice in Mathematica?

The problem is basically that my initial guess $(\theta_a = \pi/4)$ in FindRoot ends up giving me an angle that satisfies the first equation, but does not satisfy the second equation.

Answer 1

Rephrase the problem as a nonlinear optimization problem

Since $N$ is the name of a built-in function in Mathematica, we use $M$. We do the simplest nontrivial case of $M=4$. Since $A$ and $B$ are just constants, I'll set them both to 1. I'll use DiscreteDelta to take care of the $n_a^x$ and $n_a^z$ variables. We will "index" $\theta$ in the Mathematica way, as θ[a], setting θ[1]=θ[M]=0.

M=4;
A=1;
B=1;
nx[a_, M_]:=-1/Sqrt[2] DiscreteDelta[M/2-a];
nz[a_, M_]:=(1-DiscreteDelta[M/2-a]) + -1/Sqrt[2] DiscreteDelta[M/2-a];
θ[1]=0;
θ[M]=0;

Let's write the objective function to match the question but use Replace to make the $\theta_a$ periodic with $\theta_{M+1}=\theta_1$. The objective function is the square of the original function. The square has a minimum of zero where the original function is zero.

objf[a_, M_, A_, B_]:=(A(Sin[θ[a]-θ[a+1]] + Sin[θ[a]-θ[a-1]]) + 
    B(nz[a, M] Sin[θ[a]] - nx[a, M]Cos[θ[a]]))^2 /.θ[x_]->θ[(Mod[x, M]+1)]

Actually, we have $M$ equations to solve, so we need to minimize $M$ quadratics, or, equivalently, minimize their sum.

f = Sum[objf[a, M, A, B], {a, 1, M}];

We repeat the above for the constraints, except we make a list instead of a sum.

constr[a_, M_, A_, B_]:=A(Cos[θ[a]-θ[a+1]] + Cos[θ[a]-θ[a-1]]) + 
       B(nx[a, M] Sin[θ[a]] + nz[a, M]Cos[θ[a]])>0 /.θ[x_]->θ[(Mod[x, M]+1)]

c = Table[constr[a, M, A, B], {a, 1, M}];

Since $\theta_1=\theta_M=0$, our unknowns are $\theta_2, \theta_3, \theta_4, \ldots, \theta_{M-1}$. (Of course, for $M=4$, that's just $\theta_2$ and $\theta_3$.) We put the objective function, constratins, and variables into NMinimize and cross our fingers.

v=Table[θ[a], {a, 2, M-1}];
NMinimize[{f, c}, v]

 {0.255348, {θ[2] -> 4.63221*10^-23, θ[3] -> -0.261482}}

A quick Plot3D confirms this is at least reasonable.