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.

  • $\begingroup$ 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)? $\endgroup$ Dec 13, 2018 at 18:26
  • $\begingroup$ $A$ Should multiply both of the first terms. Yes, I'm using periodic boundary conditions $\endgroup$
    – MOOSE
    Dec 13, 2018 at 18:28
  • $\begingroup$ 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. $\endgroup$ Dec 13, 2018 at 18:50
  • $\begingroup$ 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/…. $\endgroup$ Dec 13, 2018 at 19:12
  • $\begingroup$ @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. $\endgroup$
    – MOOSE
    Dec 14, 2018 at 10:21

1 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.

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];

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.