0
|
14694 Views
|
19 Replies
|
6 Total Likes
View groups...
Share
GROUPS:

# How to solve a system of Non-linear equations

Posted 11 years ago
 Hi!I am trying to solve two non-linear equations, for a Tomlinson model, as follows: quauni2dx =   Quiet[(NSolve[        Subscript[V,              z]*((4*Pi*Cos[(2*Pi*y)/(Sqrt[3]*a)]*Sin[(2*Pi*x)/a])/a) +            k*(x + 0) == 0 &&                           Subscript[V,             z][(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 - Subscript[y, 0]], {x, y}, Reals] & ) /. {a -> 246/100, k -> 25/16,     Subscript[V, z] -> 5/10}]But, the Nsolve do not converge.
19 Replies
Sort By:
Posted 11 years ago
 Now, I have made some changes introducing Block, and running x -> as i, and y -> as j,, see below: $HistoryLength = 0; ClearSystemCache[]; ClearAll["Global*"]; SetDirectory["/home/marcelo/Documentos/Tomlinson/2D"]; 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} /. Take[Flatten[ 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]]; Do[ If[j == 0, c = {i, 0}]; c = Emin[i, j] ; FmolaX = -16*k*(First[c] - (i))*(10^-10); FmolaY = -16*k*(Last[c] - (j))*(10^-10); {FmolaX, FmolaY} >>> "DataFxFyTomlinsonV19.txt"; c >>> "DataMinimTomlinsonV19.txt", (*Print[ "the minima of the function related to step Y",i," X",j,":", c],*) {i, 0.0, 10.0, 0.01}, {j, 0.0, 10.0, 0.01}] ];But , an error message appear, and yet it evaluates and print the results, why? FindMinimum::lsbrak: Unable to bracket a minimum along the direction {0.,0.} from the point {3.86245,4.31621}. >> FindMinimum::lsbrak: Unable to bracket a minimum along the direction {0.,0.} from the point {3.86245,4.31621}. >> FindMinimum::lsbrak: Unable to bracket a minimum along the direction {0.,0.} from the point {3.86245,4.31621}. >> General::stop: Further output of FindMinimum::lsbrak will be suppressed during this calculation. >> FindMinimum::fdssnv: Search specification Method->PrincipalAxis without variables should be a list with 1 to 4 elements. >>ReplaceAll::reps: {FindMinimum[{{x,First[c]},{y,Last[c]}},Method->PrincipalAxis]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>ReplaceAll::reps: {FindMinimum[{{x,First[c]},{y,Last[c]}},Method->PrincipalAxis]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>FindMinimum::fdssnv: Search specification Method->PrincipalAxis without variables should be a list with 1 to 4 elements. >>ReplaceAll::reps: {FindMinimum[{{x,First[c]},{y,Last[c]}},Method->PrincipalAxis]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation. >>FindMinimum::fdssnv: Search specification Method->PrincipalAxis without variables should be a list with 1 to 4 elements. >>General::stop: Further output of FindMinimum::fdssnv will be suppressed during this calculation. >>FindMinimum::srect: Value {x,y} in search specification {x,First[c]} is not a number or array of numbers. >>FindMinimum::srect: Value {x,y} in search specification {x,First[c]} is not a number or array of numbers. >>FindMinimum::srect: Value {x,y} in search specification {x,First[c]} is not a number or array of numbers. >>General::stop: Further output of FindMinimum::srect will be suppressed during this calculation. >> Posted 11 years ago  Dear Marcelo De Cicco, could you please take a few minutes to read this tutorial about correct posting  especially of Mathematica code:How to type up a post: editor tutorial & general tipsThank you. Posted 11 years ago  Hi Moderation,I did the above instructions, (I was using, Open Suse 13.1, kernel version 13.10. , and firefox late version), but every time I published my reply, my text and codes became a terrible mess, I tried several times, and finally I had to use crt C crt V, to post my comments.Probably a bug? Posted 11 years ago  Hi Yang,After studying your graphic, I did the following code that works, or at least output the list of pairs (x,y):$HistoryLength=0;ClearSystemCache[];ClearAll["Global*"];SetDirectory["/home/marcelo/Documentos/Tomlinson/2D"];Module[{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}/.Take[Flatten[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},{y,Last}},Method->"PrincipalAxis"],7]],-2]];                                       Do[                                              If[j==0,c={0,i},c={First,Last}];                                            c=Emin[j,i]  ;                                                                          FmolaX=-16*k*(First-(j))*(10^-10);                                             FmolaY=-16*k*(Last-(i))*(10^-10);                                             {FmolaX,FmolaY}>>>"DataFxFyTomlinsonV15.txt";                                              c>>>"DataMinimTomlinsonV15.txt",                                                                                      (*Print["the minima of the function related to step  Y",i," X",j,":",c],*)  {i,0.0,10.0,0.01},                                                   {j,0.0,10.0,0.01}]                                       ];I explain the moto:the equation :-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)^2has to be minimized, which results  {x,y}, I need them as input for the next step of my "Do", they are insert as g0 , h0, respectively: Emin[g0,h0] function.There is a detail to be respect: as soon as the loop reach (i, "10")the next step (i,"0") do not should get the previous x,y, so I use the "If "command to prevent this.See that I am runnig across i: 0 to 10, thgrouh J: 0 to 10, at each step ofI I run J :0 to 10, ok?Important:   "i" corespond to h0 and "J " corresppond to g0     ( x0 is g0 and h0 is  y0, for each step -(x-x0)² and (y-yo)²)But when I graphic the results (as example the DataFxFyTomlinsonV15.txt ,results) after doing some mining dataUnprotect;ListDensityPlot[{elemFx}, Mesh -> None, InterpolationOrder -> 2, ColorFunction -> GrayLevel]Unprotect;The picture is not correct to the pic that I have to reproduce....The core is I need always the very next point (x,y) minimum, for the following step minimizing.Are my codes correct?
Posted 11 years ago
 @Cicco, Since these equations are rather complicated, I would use some graphical method to handle them first before I do the computation: I can quickly check the results, typically the distribution of the root from the following animation: Lets take a closer look at the small region around intersection Now I have some intuition about where the root might be (actually root loci ). If we are just looking for one possible result, it should be straightforward. For accuracy and efficiency, I would recommend using the optimization method. Because we are looking for the instances that both expressions equal zero, it is equivalent to find the minimum value of the square sum of  the pair of expressions, which ideally should be zero (aka least square = 0). Thus the corresponding stationary point would be one option of the roots. Particularly, FindMinimum function is the one that I would use.  To have the optimization function return a precise result, it is usually important to set a good initial point. The animation tells us that the diagnal line is not a bad choice. Also as the parameter i moves from 0 to ten, one intersection tends to move from 0 to -10. Thus I can track the root with this piece of code: (Chop clears the residues) The animation of contours of the cost functions is illustrated below (the white dot is the root/ the darker the region is, the smaller the cost function/square sum is )Here offset means how much the least squre sum is away from zero. I also heuristically adjust the initial guess to  {-i+0.1, -i+0.1} so I can have a better approximation. If you just want the numeric result, you simply remove the Manipulate wrapper and use a table function to compute "res" over a list.  Copiable code:Manipulate[ ContourPlot[{ 100/123 \[Pi] Cos[(100 \[Pi] y)/(123 Sqrt[3])] Sin[(100 \[Pi] x)/123]+(25 (x+i))/16,  1/2 ((200 \[Pi] Cos[(100 \[Pi] x)/123] Sin[(100 \[Pi] y)/(123 Sqrt[3])])/(123 Sqrt[3])+(200 \[Pi] Sin[(200 \[Pi] y)/(123 Sqrt[3])])/(123 Sqrt[3]))+(25 (y+i))/16},{x,-i-3,-i+3},{y,-i-3,-i+3}],{i,0,10,0.02}] Manipulate[  res=Chop[FindMinimum[Total[pairs[i]^2],{{x,-i+0.1},{y,-i+0.1}}]];   coord=res[[2]][[All,2]];   ContourPlot[Total[pairs[i]^2],{x,-i-2,i+4},{y,-i-2,i+4},     PlotLegends->Automatic,Epilog-> {White,PointSize[Large],Point[coord]},     PlotLabel-> Row[{"Root {x,y}=", coord, ", Offset = ", res[[1]]}],      Contours->{1,2,5,10,15,20,50,100,500,100,200,500,1000,5000},ContourStyle->None],{i,0,10,0.02},   Initialization:> (     pairs[i_]:={100/123 \[Pi] Cos[(100 \[Pi] y)/(123 Sqrt[3])] Sin[(100 \[Pi] x)/123]+(25 (x+i))/16,1/2 ((200 \[Pi] Cos[(100 \[Pi] x)/123] Sin[(100 \[Pi] y)/(123 Sqrt[3])])/(123 Sqrt[3])+(200 \[Pi] Sin[(200 \[Pi] y)/(123*Sqrt[3])])/(123 Sqrt[3]))+(25 (y+i))/16    };   ),   SynchronousUpdating->False]
Posted 11 years ago
 Shenghui,Awesome work ! I ll have take couple days to "digest" all this
Posted 11 years ago
 (*A number of changes to your code*) In[1]:= a = 1; k = 1; V = 1; c = {0.1, 0.1}; Do[c = {x, y} /. Take[Flatten[FindMinimum[-V (2 Cos[(2 Pi x)/ a] Cos[(2 Pi y)/(a Sqrt[3])] +     Cos[(4 Pi y)/(a Sqrt[3])]) + k (x - j)^2 + k (y - i)^2, {{x, First[c]}, {y, Last[c]}}]], -2];  Print[c];   Print["the minima of the fuction related to step XY ", i, j " : ", First[c]],  {i, 0, 2}, {j, 0, 2} ]Out[5]= {1.71831*^-10, 1.88021*^-10}the minima of the fuction related to step XY 001.718*^-10{0.02480, -5.15423*^-20}the minima of the fuction related to step XY 0 : 0.02480{0.05021, -1.45609*^-27}the minima of the fuction related to step XY 02  : 0.05021{1.22539*^-10, 0.02480}the minima of the fuction related to step XY 101.225*^-10{0.02490, 0.02490}the minima of the fuction related to step XY 1 : 0.02490{0.05043, 0.02521}the minima of the fuction related to step XY 12  : 0.05043{6.57703*^-10, 0.05021}the minima of the fuction related to step XY 206.577*^-10{0.02521, 0.05043}the minima of the fuction related to step XY 2 : 0.02521{0.05110, 0.05110}the minima of the fuction related to step XY 22  : 0.05110
Posted 11 years ago
 Thanks again Bill!I wil make some tests , I ll return soon.
Posted 11 years ago
 Hi again,I have done some changes in my equations, I am using FindMinimum, as I need local minima. Now I have the following issue: I need do solve x, for every Y running 0 up to10, for example. And the pair of the solutions will be the next guess to my function findminimun.Below the code I have done: InputForm[c = {0.1, 0.1}];     Do[c =          InputForm[                                                    \                                                       {x, y} /.                                                 Take[                                                            Flatten[                                                                   FindMinimum[-V (2 Cos[(2 \[Pi] x)/a] Cos[(2 \[Pi] y)/(             a Sqrt[3])] + Cos[(4 \[Pi] y)/(a Sqrt[3])]) +         k (x - j)^2 + k (y - i)^2,                                                                      \                            {{x, First[c]}, {y, Last[c]}}                                                                      \                 ]                                                                      \      ]                                                   , -2]                                                                                                                                 ]; Print[c],        Print["the minima of the fuction related to step XY ",i,j " \: ",First[c]]              {i, 0, 2}, {j, 0, 2}]  But, when it evaluates, it has a error message: {-1.9621919242424415*^-11, 9.124541851054263*^-10}  FindMinimum::nrnum: The function value {0.0625,0.0625} is not a real number at {x,y} = {{-1.96219*10^-11,9.12454*10^-10},{-1.96219*10^-11,9.12454*10^-10}}. >>  FindMinimum::nrnum: The function value {0.0625,0.0625} is not a real number at {x,y} = {{-1.96219*10^-11,9.12454*10^-10},{-1.96219*10^-11,9.12454*10^-10}}. >>  FindMinimum::nrnum: The function value {0.0625,0.0625} is not a real number at {x,y} = {{-1.96219*10^-11,9.12454*10^-10},{-1.96219*10^-11,9.12454*10^-10}}. >>  General::stop: Further output of FindMinimum::nrnum will be suppressed during this calculation. >>ReplaceAll::reps: {FindMinimum[-V (2 Cos[Times[<<4>>]] Cos[Times[<<4>>]]+Cos[4 \[Pi] y Power[<<2>>]])+k (x-j)^2+k (y-i)^2,{{x,First[c]},{y,Last[c]}}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>General::stop: Further output of$RecursionLimit::reclim will be suppressed during this calculation. >> .......
Posted 11 years ago
 Bill,I suppose  that, as I saw your code. In fact, these two expressions, are originally form a gradient (2D )apllied over an energy equation.The first result expression is concern to the i-componente and the second one to the j-component.
Posted 11 years ago
 Your two expressions each might be positive or negative. By squaring each, adding them and then searching for some minimum value this will quickly let me find a location where both expressions are very close to zero, whether each is positive or negative.
Posted 11 years ago
 Bill,I notice you square each equation, is there any reason to do that?
Posted 11 years ago
 a = 246/100; k = 25/16; V = 5/10;Table[{i, NMinimize[((V*(4*Pi*Sin[(2*Pi*x)/a]*Cos[(2*Pi*y)/(Sqrt[3]*a)]))/a + k*(x + i))^2 +  (V*((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 + i))^2,   {x, y}, MaxIterations -> 1000, Method -> "RandomSearch"]}, {i, 0, 10, 0.02}]That rapidly provides some points very close to solutions of your original system. Perhaps this will be good enough for you. If not then perhaps you can use those as starting locations to then look for more precise solutions, but if you inspect all the results you can see some values of i find very large minima.
Posted 11 years ago
 Bill,Great!I will try and then I ll send you a feedback.
Posted 11 years ago
 Well, I made clean my code, as bellow:eq = {(V*(4*Pi*Sin[(2*Pi*x)/a]*Cos[(2*Pi*y)/(Sqrt[3]*a)]))/a +     k*(x + #),        V*((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 + #)}; a = 246/100; k = 25/16; V = 5/10; I expalin my goal, I need to run the solution of this coupled system  in step of Table[i, {i, 0, 10, 0.02}, as I need the pairs x,y, for each poitn of the list, so I did the following:Quiet@NSolve[Thread[eq == 0], {x, y},   Reals & /@ Table[i, {i, 0, 10, 0.02}]]But it not converge, i think
Posted 11 years ago
 Thanks. I will try again.
Posted 11 years ago
 Move your /. inside your NSolve, provide a value for y0, add an ==0 or something similar to the end of your second equation, fix several [] that possibly should be (), discard the unused & and see if there are any more errors. Perhaps something like this:quauni2dx = Quiet[NSolve[{    Subscript[V, z]*((4*Pi*Cos[(2*Pi*y)/(Sqrt[3]*a)]*Sin[(2*Pi*x)/a])/a) + k*(x + 0) == 0,     Subscript[V, z]*((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 - Subscript[y, 0]) == 0    } /. {a -> 2.46, k -> 1.5625, Subscript[V, z] -> .5, Subscript[y, 0] -> .25}, {x, y}, Reals]]
Posted 11 years ago
 Did you try FindRoot ?Narasimham
Posted 11 years ago
 Yes, but do not converge.