Message Boards Message Boards


Using NSolve does not work : non linear equations

Posted 9 years ago
13 Replies
1 Total Likes

I have the follwing equation that I want solve:
 Block[{x, y, X = 0, Y = 0, a = 2.46, k = 1.5625},
   NSolve[{( (
          4 \[Pi] Cos[(2 \[Pi] y)/(Sqrt[3] a)] Sin[(2 \[Pi] x)/a] )/
          a) + k  (X - x)}^2 + {((
           4 \[Pi] Cos[(2 \[Pi] x)/a] Sin[(2 \[Pi] y)/(Sqrt[3] a)])/(

          Sqrt[3] a) + (4 \[Pi] Sin[(4 \[Pi] y)/(Sqrt[3] a)])/(

          Sqrt[3] a) ) + k (Y - y)}^2 == 0, {x, y}, Reals ]
 But, when it isevaluated nothing happens, just repeats the lines above:
NSolve[{(-1.5625 x +

      5.10828 Cos[1.47463 y] Sin[2.55414 x])^2 + (-1.5625 y +

      2.94927 Cos[2.55414 x] Sin[1.47463 y] +

      2.94927 Sin[2.94927 y])^2} == 0, {x, y}, Reals]
Is it possible that couple equations does not have a set of solutions?
POSTED BY: Marcelo De Cicco
13 Replies
Posted 9 years ago
Hi Marcelo,

Thanks for the information. I have some questins and comments:

Here is what I get for your potential function, with all constants = 1, and (X,Y) = (0,0).

First, in looking at your potential, I see that there are two terms:
      The first is periodic in x and y with a hexagonal cell structure having a period defined by a lattice constant a.

      The second is parabolic upward, with curvature k and center at (X,Y).

I note that the constant V multiplies only the periodic term, which is perhaps not as you intend since it produces inconsistent units.

With the periodic term negated as shown, the result is that the periodic term has a minima at (x,y) = (0,0), ant of course many more.

The parabolic term with (X,Y) = (0,0) has a single minimum at (x,y) = (0,0) also.

The iteration you describe seems to be interating the position of the parabolic term with the goal of minimizing the minimum of the resulting potential. Two comments on that:
  1. Finding the minimum of the combined potential, with an off-center parabolic term, is a global minimization over a function with a great many local minima. This can be very troublesome.
  2. If indeed only the periodic term is negated as shown, the minimization problem has a solution by inspection which occurs whenever the parabolic term is centered on a cell. For example  (X,Y) = (0,0) is a solution.
So, generally speaking, I do not really understand the purpose of that algorith. (And indeed, I might not understand.)

I would like to suggest another possible approach to fitting your theoretical potential to the experimental data, and recovering matching plots.
  1. I assume you have Fx(x,y) and Fy(x,y) as experimental data, and can assure their coordinate systems coincide.Create an Interpolating function for each.
  2. From this you can construct an equivalent experimental potential function by path integration. Choose any convenient point in the field as the zero of potential and the potential at any other point is given by the path integral of  -F dl. This can be easily done for any point by choosing a path first parallel one axis, and then the other. For example, Integral -F dl = Integral -Fx dx + Integral -Fy dy along the two legs of the path. This could be done in a Table using Nintegrate.
  3. This experimenatal data list (appropriately flattened) can be used as data to which the potential function is curve fit.
  4. Before trying to fit the data, look at the periodic nature of the data and calculate a good guess for the lattice constant a. Also the best guess you can for the other parameters. You will also need to add a constant fit parameter to V (V+c) to the potential model to allow for matching the choice of a zero of potential.
  5. Now use NonlinearModelFit to fit your theoretical potential function to the potential table you derived from the force data. If you have chosen close starting points this should work well.
  6. To get to the pictures, generate Fx(x,y) = -d/dx V(x,y) and Fy(x,y) = -d/dy V(x,y) where the derivatives are partial and V(x,y) is your fit model.
Kind regards,
POSTED BY: David Keith
Posted 9 years ago
Thank you, Marcelo.  I will think about this.
Kind regards,
POSTED BY: David Keith

I will be more specific, starting the problem since begining:

1)I have these images (map of force due a potential, Fx and Fy) that I need to replicate the last two , using mathematica ( my  adviser wanted to use Fortran..., but  I am fighting to consolidate Mathematica coding in my department )


The potential that I need to Minimize at each step is:
-V *(2 Cos[(2 Pi x)/a]*

  Cos[(2 Pi y)/(a Sqrt[3])] +  Cos[(4 Pi y)/(a Sqrt[3])]) + (k/2)*(X- x)^2 + (k/2)*(Y - y)^2
So, when I get the minima, or (x,y) positions, they will be used to the next step. I  use them like that:
 first step:
  i(=X) = 0 and j (=Y)=0,  -> potential minimized -> x,y calculated.
 second step:
 i= 0 j=1, potential minimized using the last (x,y).....-> "new" x,y calculated..
 and so on...
And I am going collecting at each step the pairs (x,y) , that will be used to calculate the Map Forces Fx and Fy.

At the final I have an output.txt with all x,y, collected. and I use x for calculate Fx and y for calculate Fy.

and finally I make a table to plot Fx x X(or i) and Fy x Y(or j)

obs: "F = k* x" (oscilator)

well, and  I did a lot of tests (since august....) and the best images I could get (using contourplot) were (for Fx, and Fy, respectively)

The parameters:
(*a= 2.46 \[Angstrom]  = 2.46 10^-10m*)

(*k= 25 N/m  as 1 joule = N.m, so , N/m= eV/(16 \[Angstrom]^2)*)

And the code  is that:
 $HistoryLength = 0;
 Block[{x, y, i, j, FmolaX, FmolaY, c = {0, 0}, a = 2.46, k = 1.5625,

   V = 0.5},

       Emin[go_, ho_] :=

   Module[{g = go, h = ho}, {x, y} /.


       N[FindMinimum[-V *(2 Cos[(2 Pi x)/a]*

              Cos[(2 Pi y)/(a Sqrt[3])] +

             Cos[(4 Pi y)/(a Sqrt[3])]) + (k/2)*(x - g)^2 + (k/

             2)*(y - h)^2, {{x, First[c]}, {y, Last[c]}},

         Method -> "PrincipalAxis"], 7]], -2]];


                                               If[j == 0, c = {0, i}];         (*<- to avoid contour problems in the beginning of the step, I can´t use the last x,y results, for the next i*)                              

   c = Emin[j, i]  ;                                                                           

   FmolaX = -16*k*(First[c] - (j))*(10^-10);                                              

   FmolaY = -16*k*(Last[c] - (i))*(10^-10);
     {FmolaX, FmolaY} >>>  "DataFxFyTomlinsonV20.txt";

                    c >>> "DataMinimTomlinsonV20.txt",


   "the minima of the function related to step  Y",i," X",j,":",

   c],*)  {i, 0.0, 10.0, 0.1},   {j, 0.0, 10.0, 0.1}]

If you can  see, I did steps of 0.1, and each step I used the last pair of {x,y},of the minimized potential.

As my images were not so good, mainly Fy, I tried aplying the gradient :
F=-\[Del](x,y) Subscript[E, tot]

So , I got two eqautions from the potential, and then I tested NSolve for Fx=0 and Fy = 0.

Details, When I tested the model 1d worked fine! See below:
 equauni =
   Quiet[NSolve[(2 \[Pi] )/ a Subscript[V, z]  Sin[(2 \[Pi] x)/a] +
          k  (x - #) == 0, x, Reals] & /@
      Table[i, {i, 0, 10, 0.02}] /. {a -> 246/100, k -> 25/16,
      Subscript[V, z] -> 5/10}] >> datateste1

equauni2 = Flatten[Take[{x} /. equauni // Simplify, All, 1]];

dataequani = Transpose[{Table[i, {i, 0, 10, 0.02}], equauni2}];

Map[Length, {equauni2}];

Map[Length, {Table[i, {i, 0, 10, 0.02}]}];

ListPlot[Tooltip[Transpose[{Table[i, {i, 0, 10, 0.02}], equauni2}]],

PlotStyle -> Green,

AxesLabel -> {Subscript[x, 0][\[Angstrom]], x[\[Angstrom]]}]
And calculating the Force:
    d_] := -k  (x -
      d)*(10^-10);(*converted to  Newton (N)*)
   Transpose[{  10^-10 equauni2,

    Fmola[equauni2, Table[i, {i, 0, 10, 0.02}]] /. {a -> 246/100,

      k -> 25/1, Subscript[V, z] -> 5/10}}]],

AxesLabel -> {x["m"], F["N"]}, AxesOrigin -> {0, 0},

PlotRange -> All, ImageSize -> Large, PlotStyle -> Green, Mesh -> All]
I hope, this explanation can be more useful.
POSTED BY: Marcelo De Cicco
Posted 9 years ago
Hello Marcelo and Bruce,

I understand what you mean, Bruce.

But I still do not really understand the problem. This nonlinear implicit transcendental equation in x and y has been separately described as two equations to be solved, and as an objective function to be minimized. Two days ago (above) I used FindMinimum on a grid of starting values to hunt the zeros of this function, assuming it is an objective function in two dimensions (R2->R1),  and found and plotted a number of them, indicating that if this is an objective function to be minmized, it has a number of local minima.

So I am really not sure what is of interest.

POSTED BY: David Keith

There are a few non-algebraic special cases, mostly things not far from algebra.
In[1]:= NSolve[E^x - x == 7, x, Reals]                                                 

Out[1]= {{x -> -6.99909}, {x -> 2.22154}}
Someday when I have the time, I will try to show sampling solutions and interpolating.
This is not a good week for that.  Maybe someone else can take it on?
POSTED BY: Bruce Miller
In, under Details and Options,
most of the way down,

    NSolve deals primarily with linear and polynomial equations.

NSolve appears to be trying to figure out the Sin and Cos terms, but I would be surprised
if you can get symbolic solutions for x and y.

It would be a little work than I have time for now, but one could use FindRoot to find the y that solves
the equation for a given x, do that for many xs,  combine those xs and ys into a list {.., {xi,yi}, ..}, then Interpolate that.
POSTED BY: Bruce Miller

I see your point, but as the quotes read  "primarily..." : is not excludent.

I do not understand very well your suggestion, can you be more explicit? 

POSTED BY: Marcelo De Cicco

Thanks, I will see your codes, but previously, I need  a set of solutions for x and y, as I will make steps (I walk in x and y, using always the previousx, y results) I need exactly  very next pairs x,y .

Nsolve I will give me a set of solutions
POSTED BY: Marcelo De Cicco
We can also note that the surface goes negative. We ca therefore square everything and find minima:
 NMinimize[(((4 \[Pi] Cos[(2 \[Pi] y)/(Sqrt[3] a)] Sin[(2 \[Pi] x)/a])/
        a) + k (X -
          x)^2 + (((4 \[Pi] Cos[(2 \[Pi] x)/
                a] Sin[(2 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a) + (4 \[Pi] Sin[(4 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a)) + k (Y - y))^2)^2 /. {X -> 0, Y -> 0,
    a -> SetPrecision[ 2.46, Infinity],
    k -> SetPrecision[1.5625, Infinity]}, {x, y}]
(*{2.08309*10^-30, {x -> 1.6586, y -> 0.0891708}}*)
The following figure helps to narrow down the position of the roots:
Plot3D[(1.5625` x -
    5.108280737544379` Cos[1.474633629458714` y] Sin[
      2.5541403687721895` x])^2 + (1.5625` y -
    2.949267258917428` Cos[2.5541403687721895` x] Sin[
      1.474633629458714` y] -
    2.949267258917428` Sin[2.949267258917428` y]), {x, -5, 5}, {y, -5,
   5},  PlotRange -> {All, All, {-10, 0.0}} ]

POSTED BY: Marco Thiel

I noticed that some of the parenthesis do not seem to be placed correctly. There are to sets of curly one that do not seem to work. If you correct that it is easy to find solutions:
  Simplify[((4 \[Pi] Cos[(2 \[Pi] y)/(Sqrt[3] a)] Sin[(2 \[Pi] x)/a])/
        a) + k (X -
          x)^2 + (((4 \[Pi] Cos[(2 \[Pi] x)/
                a] Sin[(2 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a) + (4 \[Pi] Sin[(4 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a)) + k (Y - y))^2 == 0] /. {X -> 0., Y -> 0.,
    a -> 2.46, k -> 1.5625}, {x, y}, Reals]
(*{{x -> 0, y -> 0}}*)

 The SolveAlways command might also help:
Simplify[((4 \[Pi] Cos[(2 \[Pi] y)/(Sqrt[3] a)] Sin[(2 \[Pi] x)/a])/
       a) + k (X -
         x)^2 + (((4 \[Pi] Cos[(2 \[Pi] x)/
               a] Sin[(2 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
              3] a) + (4 \[Pi] Sin[(4 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
              3] a)) + k (Y - y))^2 == 0] /. {X -> 0., Y -> 0.,
   a -> 2.46, k -> 1.5625}, y]

This gives a rather long output. 

This command
ContourPlot[(1.5625` x -
    5.108280737544379` Cos[1.474633629458714` y] Sin[
      2.5541403687721895` x])^2 + (1.5625` y -
    2.949267258917428` Cos[2.5541403687721895` x] Sin[
      1.474633629458714` y] -
    2.949267258917428` Sin[2.949267258917428` y])^2, {x, -5,
  5}, {y, -5, 5}, PlotLegends -> Automatic, Contours -> 50]

Might help to narrow down the search.

POSTED BY: Marco Thiel
Posted 9 years ago
Hi Marcelo,

Two issues first: 
1) The curly bracekets are used to denote lists. Perhaps where they appear in the equation you should be using ordinary parentheses as outer brackets?
2) If the curly brackets are so replaced, the result is a single equation in two variables, so the problem is under constrained. Did you mean there to be two equations here? Or are you trying to solve the implicit equation in x and y for some explicit y of x? If so, Solve is what is needed.

POSTED BY: David Keith
Hi David,

1) Sometimes It can be used curly brackets in Nsolve.

2) And yes, they are two equations, derived from:
-V (2 cos((2 \[Pi] x)/a) cos((2 \[Pi] y)/(Sqrt[3] a))+cos((4 \[Pi] y)/(Sqrt[3] a)))+1/2 k (X-x)^2+1/2 k (Y-y)^2
I need all minimized solutions form there, as it is an objective functions to be optimezed, but when I usedFind MInimum, I had only two solutions for x and y, but I need a set os solutons.
POSTED BY: Marcelo De Cicco
Posted 9 years ago
Hi Marcelo,

I still do not see two equations. There is only one equal sign, and it is not between two conforming lists, or a pair of complex numbers. I do see that the syntax is allowed. Mathematica appears to treat {a+b}^2 = { (a+b)^2 }. So I used the equation given as the left hand side being a merit function to be minimized. Ideally to find its zeros.

Following Marcos, I note that the central region is of interest. The code below defines the merit function and values as I think you have given them. It then uses FindMinimum with a grid of starting points to generate a hopefully complete list of the numerical zeros. Chop is used to make very small numbers zero. Union eliminates duplicates and provides a list of the minimum coordinates, which are plotted. Note that Mathematica's NMinimize will also look for a global mimimum, and return one of those found with the grid search. Note that with the constant values you've given, {0,0} is always a zero.

Am I seeing this correctly?
 (* the merit function *)
 mf = Simplify[{((4 \[Pi] Cos[(2 \[Pi] y)/(Sqrt[3] a)] Sin[(2 \[Pi] x)/
               a])/a) +
         k (X - x)}^2 + {((4 \[Pi] Cos[(2 \[Pi] x)/
                a] Sin[(2 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a) + (4 \[Pi] Sin[(4 \[Pi] y)/(Sqrt[3] a)])/(Sqrt[
               3] a)) + k (Y - y)}^2][[1]];
 (* constants *)values = {X -> 0., Y -> 0., a -> 2.46, k -> 1.5625};

(* for the values given for the constants there is always a zero at \
the origin *)mf /. x -> 0 /. y -> 0;

(* function to minimize the mf, starting at {x0,y0} *)
minimizeMf[x0_, y0_] := FindMinimum[mf /. values, {{x, x0}, {y, y0}}]

(* build a table of local minimizations on a grid *)
tries = Flatten[
    Table[minimizeMf[x0, y0], {x0, -5, 5, .5}, {y0, -5, 5, .5}], 1] //

(* select those with zero value, eliminate duplicates, and keep only \
the {x,y} values *)
zeros = Union[{x, y} /. Select[tries, #[[1]] === 0 &][[All, 2]],
   SameTest -> (Norm[#1 - #2] < 10^-10 &)];

(* plot the found zeros on a contour plot *)
ContourPlot[mf /. values, {x, -5, 5}, {y, -5, 5}, Contours -> 50,
PlotLegends -> Automatic, Epilog -> {Red, Point /@ zeros}]

In[8]:= (* Mathmaticas global minimization would find one of these *)

NMinimize[{mf /. values, -5 < x < 5 && -5 < y < 5}, {x, y}]

Out[8]= {1.75475*10^-18, {x -> 0.548069, y -> 0.949283}}
POSTED BY: David Keith
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract