# Finding the minimum points of a discrete trigonometric equation

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.

• Should $A$ multiply both of the first two terms in the first equation? Also, is $\theta_{N+1}$ the same thing as $\theta_1$ (i.e., periodic boundary conditions)? – Michael Seifert Dec 13 '18 at 18:26
• $A$ Should multiply both of the first terms. Yes, I'm using periodic boundary conditions – MOOSE Dec 13 '18 at 18:28
• You could always generate a large number of roots of the first system using randomized initial guesses and then select all those for which the second conditions hold. It's a bit kludgy, but it might work. – Michael Seifert Dec 13 '18 at 18:50
• This may be a naive suggestion, but can you rephrase the problem as a constrained optimization problem? E.g. minimize the square of the LHS subject to the given constraints. If so, check out Mathematica's overview of Constrained Optimization: reference.wolfram.com/language/tutorial/…. – Robert Jacobson Dec 13 '18 at 19:12
• @RobertJacobson it might be! Do you know if it is easy to minimize this problem using for instance NMinimize? I'm having some difficulty due to that the equations are coupled. – MOOSE Dec 14 '18 at 10:21

# 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 θ=θ[M]=0.

M=4;
A=1;
B=1;
nx[a_, M_]:=-1/Sqrt DiscreteDelta[M/2-a];
nz[a_, M_]:=(1-DiscreteDelta[M/2-a]) + -1/Sqrt DiscreteDelta[M/2-a];
θ=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, {θ -> 4.63221*10^-23, θ -> -0.261482}}


A quick Plot3D confirms this is at least reasonable. 