Equations and Systems of Equations

Linear and nonlinear equations

There are two groups of equations and systems of equations: linear and nonlinear. This division holds for both equations for unknown quantities (numbers) that are considered in this chapter and equations for unknown functions such as differential equations that will be considered later. Linear equations contain only zero and first combined powers of unknowns, such as the equation

ax==b

and the system of equations

ax+by==c
dx+ey==f.

Here x and y are unknowns and other symbols are parameters or numbers. Equations containing higher combined powers of unknowns such as (x^2), xy, or functions of unknowns such as Sin[x], are nonlinear.

Linear equations and systems of equations are analytically solvable and the solution, if it exists, is unique. Nonlinear equations can have several solutions and analytical solutions are available in a limited number of cases such as quadratic, cubic, and quartic equations

(ax^2)+bx+c==0
(ax^3)+(bx^2)+cx+d==0
(ax^4)+(bx^3)+(cx^2)+dx+e==0,

as well as some irrational equations

(Sqrt[x-a])+(Sqrt[x-b])==c,

some trigonometric equations

Sin[x]==Cos[2x]

and other equations. All these equations are algebraic because they can be simplified to polynomial equations

If equations contain different functions of unknowns at the same time, such as

x==Cos[x]

they are called transcedental equations because they are non-algebraic and never can be solved analytically. In many cases solution of nonlinear equations must be done numerically and special care should be taken in the case of many solutions.

Solution of equations with Mathematica

Linear equations

Both analytical and numerical solutions of equations are output by Mathematica in the form of rules

Solve[a x+b==0,x]

that with the help of replacement can be further used to find the resulting expression for x

x/.Solve[a x+b==0,x]

or better

x/.Solve[a x+b==0,x][[1]]

or to define x as a function of the parameters

xab[a_,b_]=x/.Solve[a x+b==0,x][[1]];
xab[a,b]

Nonlinear algebraic equations

Note that the solution of an equation is a list of rules corresponding to all its roots. Quadratic equation has two roots,

sol=Solve[a (x^2)+b x+c==0,x]

so that sol is a list of two elements. The roots of this equations form the list

x/.sol

and the individual roots can be addressed as

x/.sol[[1]]
x/.sol[[2]]

The following irrational equation can be reduced to a quadratic equation by squaring it two times. The resulting quadratic equation has two roots, only one of which is the solution of the initial equation, why the other solution is parasite. This is why this equation is still algebraic.

sol=Solve[(Sqrt[x-a])+(Sqrt[x-b])==c,x]

The simplified expression for this solution is

sol={{x->((Simplify[(a^2)-2 a b+(b^2)]+Simplify[2 a (c^2)+2 b (c^2)+(c^4)])/(4 (c^2)))}}

Below is a particular case of the above equation. Since x is the only symbol here, it is the unknown and it does not have to be indicated in the Solve command

Solve[(Sqrt[2x])+(Sqrt[x+1])==(Sqrt[2])]
N[%]

Some equations can be solved analytically in terms of complex numbers

Solve[(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1])==3]

This solution is real but Mathematica cannot simplify this expression to the explicitly real form. Conversion to the numerical form shows that the solution is real indeed, up to the imaginary numerical noise term

N[%]

The latter can be eliminated with the help of the Chop command

Chop[%]

In such cases it is advisable to illustrate the equation graphically

Plot[{(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1]),3},{x,1,1.5},PlotRange->All]

Indeed, the solution above is OK. More practical is to force Mathematica to solve this equation numericaly from the beginning

Solve[(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1])==3.]

For more complicated equations there can be an output saying that Mathematica has reduced the equation to a polynomial equation

Solve[(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1])+(Sqrt[x+2])==5]

that can be solved only numerically

N[%]

Again, numerical solution can be done from the beginning

Solve[(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1])+(Sqrt[x+2])==5.]

In this case one can use the NSolve command that works on polynomial equations, according to Mathematica's help

NSolve[(Sqrt[x-1])+(Sqrt[x])+(Sqrt[x+1])+(Sqrt[x+2])==5]

Although this equation is not polynomial, Mathematica can reduce it to a polynomial equation, this is why NSolve works here.

Many trigonometric equations can be solved by their reduction to polynomial equations

Solve[Sin[x]==Cos[2x],x]
N[%]

Indeed, this equation is equivalent to a quadratic equation for Sin[x] because

Cos[2x]=((Cos[x])^2)-((Sin[x])^2)=1-2((Sin[x])^2)

This formula can be obtained with Mathematica by

TrigExpand[Cos[2x]]/.((Cos[x])^2)->1-((Sin[x])^2)

The solution of the resulting quadratic equation with respect to Sin[x] has the form

Solve[Sin[x]==1-2 ((Sin[x])^2),Sin[x]]

This output constitutes two simple triginometric equations for x that can be solved to give the result above, x = -π/2, x = π/6, and x = 5π/6. Since the original equation is algebraic, also NSolve (applicable to polynomials) does the job

NSolve[Sin[x]==Cos[2x],x]

Both Solve and NSolve do not take into account the periodicity of the trigonometric function and thus lose some solutions. Complete solutions can be obtained by Reduce

Reduce[Sin[x]==Cos[2x],x]
N[%]

Another example of an algebraic equation is

Solve[(E^x)-(E^(2x))+1==0]

One of the two solutions is real and the other is complex. This is easy to understand because this equation is a quadratic equation for (e^x). An attempt to solve this equation for (e^x)

Solve[(E^x)-((((E^x)))^2)+1==0,(E^x)]

does not work, although an equation was solved for Sin[x] above. Well, with the substitution y = (e^x) one obtains

Solve[y-(y^2)+1==0]

The first solution is negative, although exponential of a real argument is positive. This, looking for real solutions of our equations, one has to drop it.

Transcedental (non-algebraic) equations

Transcedental equations cannot be solved with Solve, NSolve, and Reduce

Solve[Cos[x]==x]

NSolve[Cos[x]==x]

Reduce[Cos[x]==x]

The only way to solve such equations is to use the essentially numerical procedure FindRoot

FindRoot[Cos[x]==x,{x,0}]

0 in this command is the starting value of x. Since this equation has a single solution

Plot[{Cos[x],x},{x,-3,3}]

the choice of the starting point is irrelevant. If the equation has several solutions, usually (but not always!) FindRoot finds the solution closest to the starting point. For instance, the equation

Cos[x]==x/5

has three solutions:

Plot[{Cos[x],x/5},{x,-7,7}]

FindRoot[Cos[x]==x/5,{x,0}]

Here FindRoot got stuck in the vicinity of the near-solution x ≈ 6.

FindRoot[Cos[x]==x/5,{x,2}]

This is the correct positive solution.

FindRoot[Cos[x]==x/5,{x,-1}]

Also a correct solution.

FindRoot[Cos[x]==x/5,{x,-5}]

But

FindRoot[Cos[x]==x/5,{x,-6}]

gives a solution that is not the closest to the starting point. One can plot the dependence of the solution on the starting point

xx0[x0_]:=x/.FindRoot[Cos[x]==x/5,{x,x0}]
Plot[xx0[x0],{x0,-5,2}]

Here one can see all three roots plus the wrong root x ≈ 6. One can eliminate the wrong root by bracketing the search interval to, say,  {-5,2} in the FindRoot command

xx0[x0_]:=x/.FindRoot[Cos[x]==x/5,{x,x0,-5,2}]
Plot[xx0[x0],{x0,-5,2}]

There are functions that make problems in the standard version of FindRoot if the starting point is not sufficiently close to the root:

f[x_]=(x/(1+(x^2)));
Plot[f[x],{x,-10,10}]

FindRoot[f[x]==0,{x,0.7}]

is the correct solution but already

FindRoot[f[x]==0,{x,0.8}]

does not work, as well as

FindRoot[f[x]==0,{x,5}]

In such situations one can use {x, (Subscript[x, 1]), (Subscript[x, 2])} in the FindRoot statement. Mathematica's help says that in this case FindRoot does not use analytical derivatives of f[x] (see below on the role of derivatives) and computes derivatives numerically using two different values of x, the initial two values being (Subscript[x, 1]) and (Subscript[x, 2]). It also appears that if (Subscript[x, 1]) and (Subscript[x, 2]) are bracketing the root, the latter is found without problems:

FindRoot[f[x]==0,{x,-100,300}]

That is a correct result up to numerical errors.

FindRoot works iteratively. The number of steps and function evaluations needed to reach the root with the required precision can be output via StepMonitor and EvaluationMonitor:

s=0;e=0;
FindRoot[Cos[x]==x/5,{x,1},StepMonitor:>s++,EvaluationMonitor:>{e++,Print[e,"   ",x]}]
Print["Steps = ",s,"    Evaluations = ",e]

or

s=0;e=0;
FindRoot[f[x]==0,{x,-2,200},StepMonitor:>s++,EvaluationMonitor:>{e++,Print[e,"   ",x]}]
Print["Steps"->s,"    Evaluations"->e]

Numerical method of solving transcedental equations

The basic method for solving transcedental equations of the type F[x] == 0 is the Newton method, in which the function is linearized around the starting point (Subscript[x, 0])

(Subscript[F, lin])[x]=F[(Subscript[x, 0])]+F'[(Subscript[x, 0])](x-(Subscript[x, 0]))

and the root of the (Subscript[F, lin])[x] is found as

(Subscript[x, 1])=(Subscript[x, 0])-((F[(Subscript[x, 0])])/(F'[(Subscript[x, 0])])).

Taking (Subscript[x, 1]) as the next starting point, one repeats this procedure iteratively until |(Subscript[x, 1])-(Subscript[x, 0])|<δ. Below is the procedure code. The first argument of the Module statement is the list of local variables. As we need xList as an output, we keep it global. Also the maximal number of iterations kMax is kept global.

FindRootNewton[Func_,xStart_]:=Module[{del,k,x0,x1,FDer},
FDer[x_]=(∂_x)Func[x];
x1=xStart;  x0=xStart+1;   k=1;        del=(10^(-16));
xList={x1};    
While[Abs[x1-x0]>del&&k<kMax,
k++;
x0=x1;
x1=x0-((Func[x0])/(FDer[x0]));
xList=Append[xList,x1]
 ];
x1];

FindRootNewton is programmed as a function, the output is the value x1, not a rule. Its application is the following:

F[y_]=y+(y^2);   kMax=30;
FindRootNewton[F,1.]

For a comparison,

FindRoot[F[x],{x,1.}]

Number of iterations and iteration results in FindRootNewton:

Length[xList]
xList

One can see that the convergence becomes very fast when x becomes close to the root.

ListLogPlot[xList]

Here is the plotting procedure visualizing how the Newton's method works. kShow is the maximal number of iterations to show.

ShowHowNewtonWorks[kShow_,xLeft_,xRight_]:=Show[
Plot[F[x],{x,xLeft,xRight},PlotRange->All,PlotStyle->{Thick},AspectRatio->1],
Graphics[{Dashed,Arrowheads[Medium],Table[Arrow[{{xList[[k]],0},{xList[[k]],F[xList[[k]]]}}],{k,1,kShow}]}],
Graphics[{Arrowheads[Medium],Table[Arrow[{{xList[[k]],F[xList[[k]]]},{xList[[k+1]],0}}],{k,1,kShow}]}],
Graphics[Disk[{First[xList],0},0.01(xRight-xLeft)]],
Graphics[{Red,Disk[{Last[xList],0},0.01(xRight-xLeft)]}]
];

Now visualize

kShow=3;
ShowHowNewtonWorks[kShow,-1.5,1]

If, for the same function, one takes the starting point on the other side of the root, the situation can become dangerous. For the starting value -0.4 the procedure still converges to the same root.

F[y_]=y+(y^2);   kMax=30;
FindRootNewton[F,-0.4]
kShow=4;
ShowHowNewtonWorks[kShow,-1.5,1]

But for the starting value -0.6 and less the procedure converges to another root.

F[y_]=y+(y^2);   kMax=30;
FindRootNewton[F,-0.6]
kShow=4;
ShowHowNewtonWorks[kShow,-2,0.5]

For the function considered above, x/(1+(x^2)) and similar, the Newton method diverges for most of starting points

F[x_]=((20x)/(1+(x^2)));   kMax=3;
FindRootNewton[F,2.]
kShow=2;
ShowHowNewtonWorks[kShow,-5,15]

or even

F[x_]=((4x)/(1+(x^2)));   kMax=3;
FindRootNewton[F,0.6]
kShow=2;
ShowHowNewtonWorks[kShow,-2,2]

But if the starting point is closer to the root, the procedure converges:

F[x_]=((4x)/(1+(x^2)));   kMax=5;
FindRootNewton[F,0.5]
kShow=3;
ShowHowNewtonWorks[kShow,-2,2]

For such "pathological" functions, one can apply at first the slower converging bisection method (if the root is bracketed) and then switch to the fast-converging Newton method closer to the root. Professional-grade nonlinear-equation solvers usually employ both methods with the internal monitoring of the convergence.

Systems of equations

Generally, for a system of equations to be solvable, the number of unknowns should be equal to the number of equations. 

If the number of equations is smaller than the number of unknowns, this system of equations is called underdetermined. In this case one should look for additional equations.

If the number of equations is larger than the number of unknowns, this system of equations is called overderdetermined. In this case one should check whether some equations follow from the others and thus can be dropped. If the situation persists, something is fundamentally wrong with the model.

Below we consider only systems of equations with equal numbers of equations and unknowns.

Systems of linear equations

System of linear equations can be solved except for the cases when equations are contradictory, such as

x+y==0
x+y==1

or

x+y+z==0
x+y+2z==1
2x+2y+3z==2

Here adding the first and second equation, one obtains 2 x + 2 y + 3 z == 2 that contradicts the third equation. In all these cases the determinant of the system of equations is zero (see below).

Consider a rigid rod of mass M and length L lying on two supports. Support 1 is located at the distance (Subscript[x, 1]) from the left end of the rod and support 2 is located at the distance (Subscript[x, 2]) from the left end of the rod. Let us find the upward reaction forces (Subscript[F, 1]) and (Subscript[F, 2]) acting on the rod from each support.

There are two linear equations

(Subscript[F, 1])+(Subscript[F, 2])-M g==0    (* Sum of all forces is zero *)
(Subscript[F, 1])(Subscript[x, 1])+(Subscript[F, 2])(Subscript[x, 2])-M g L/2==0  (* Sum of all torques (here with respect to the left end) is zero *)

This system of equation can be solved with Mathematica

Solve[{(Subscript[F, 1])+(Subscript[F, 2])-M g==0,(Subscript[F, 1])(Subscript[x, 1])+(Subscript[F, 2])(Subscript[x, 2])-M g L/2==0},{(Subscript[F, 1]),(Subscript[F, 2])}]//Simplify

or by hand. The human form of this result is

(Subscript[F, 1])=((L/2-(Subscript[x, 2]))/((Subscript[x, 1])-(Subscript[x, 2])))M g,       (Subscript[F, 2])=((L/2-(Subscript[x, 1]))/((Subscript[x, 2])-(Subscript[x, 1])))M g.

In the case of three or more supports, the number of unknowns is the number of supports, whereas there are still two equations. Thus the problem becomes underdetermined. It can be made well determined if one renounces the model of the rigid rod and takes into account its elasticity. In this case, however, the problem becomes much more complicated.

The problem with two supports is very simple and it does not require Mathematica. Let us find the effective resistance of the so-called bridge consisting of five resistors.

The currents in the nodes are related by the obvious Kirchhof's equations

I==(Subscript[I, 1])+(Subscript[I, 3])
(Subscript[I, 1])==(Subscript[I, 2])+(Subscript[I, 5])
(Subscript[I, 3])+(Subscript[I, 5])==(Subscript[I, 4])

while the equation for the rightmost node (Subscript[I, 2])+(Subscript[I, 4])==I should be discarded because it follows from the above three equations. Ohm's law for the three ways from left to right have the form

(Subscript[R, 1])(Subscript[I, 1])+(Subscript[R, 2])(Subscript[I, 2])==V
(Subscript[R, 3])(Subscript[I, 3])+(Subscript[R, 4])(Subscript[I, 4])==V
(Subscript[R, 1])(Subscript[I, 1])+(Subscript[R, 5])(Subscript[I, 5])+(Subscript[R, 4])(Subscript[I, 4])==V

One can write down more equations of this type but they follow from the above equations. For the given total voltage V, one finds 6 currents from the 6 equations and then defines the effective resistance as

R=(V/I)

Solution of this problem by hand is cumbersome. One can eliminate some currents from the first three equations but then one has to solve the remaining system of three equations with three unknowns. Mathematica solves this problem as

sol=Solve[{
II==(Subscript[II, 1])+(Subscript[II, 3]),
(Subscript[II, 1])==(Subscript[II, 2])+(Subscript[II, 5]),
(Subscript[II, 3])+(Subscript[II, 5])==(Subscript[II, 4]),
(Subscript[R, 1])(Subscript[II, 1])+(Subscript[R, 2])(Subscript[II, 2])==V,
(Subscript[R, 3])(Subscript[II, 3])+(Subscript[R, 4])(Subscript[II, 4])==V,
(Subscript[R, 1])(Subscript[II, 1])+(Subscript[R, 5])(Subscript[II, 5])+(Subscript[R, 4])(Subscript[II, 4])==V
},{II,(Subscript[II, 1]),(Subscript[II, 2]),(Subscript[II, 3]),(Subscript[II, 4]),(Subscript[II, 5])}
]

and the effective resistance is given by

Simplify[(1/(II/.sol[[1]]))V]

(V has been put after the replacement because otherwise the command does not work. Also the assignment R = ... does not work.). The human form of the result is (Excersize: Try to obtain this formula with Mathematica - could be difficult or impossible)

R= ((((Subscript[R, 1])+(Subscript[R, 2]))((Subscript[R, 3])+(Subscript[R, 4]))(Subscript[R, 5])+((Subscript[R, 1])+ (Subscript[R, 3]))(Subscript[R, 2]) (Subscript[R, 4]) +(Subscript[R, 1]) (Subscript[R, 3])((Subscript[R, 2])+(Subscript[R, 4])) )/(((Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, 3])+(Subscript[R, 4]))(Subscript[R, 5])+((Subscript[R, 1])+(Subscript[R, 3]))((Subscript[R, 2])+(Subscript[R, 4])))).

This formula is symmetric with respect to {(Subscript[R, 1])⇌(Subscript[R, 2]), (Subscript[R, 3])⇌(Subscript[R, 4])} and {(Subscript[R, 1])⇌(Subscript[R, 3]), (Subscript[R, 2])⇌(Subscript[R, 4])}. Check:

(((Subscript[R, 2]) ((Subscript[R, 4]) (Subscript[R, 5])+(Subscript[R, 3]) ((Subscript[R, 4])+(Subscript[R, 5])))+(Subscript[R, 1]) ((Subscript[R, 2]) ((Subscript[R, 3])+(Subscript[R, 4]))+(Subscript[R, 4]) (Subscript[R, 5])+(Subscript[R, 3]) ((Subscript[R, 4])+(Subscript[R, 5]))))/((Subscript[R, 3]) (Subscript[R, 4])+(Subscript[R, 3]) (Subscript[R, 5])+(Subscript[R, 4]) (Subscript[R, 5])+(Subscript[R, 2]) ((Subscript[R, 3])+(Subscript[R, 5]))+(Subscript[R, 1]) ((Subscript[R, 2])+(Subscript[R, 4])+(Subscript[R, 5]))))==((((Subscript[R, 1])+(Subscript[R, 2]))((Subscript[R, 3])+(Subscript[R, 4]))(Subscript[R, 5])+((Subscript[R, 1])+ (Subscript[R, 3]))(Subscript[R, 2]) (Subscript[R, 4]) +(Subscript[R, 1]) (Subscript[R, 3])((Subscript[R, 2])+(Subscript[R, 4])) )/(((Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, 3])+(Subscript[R, 4]))(Subscript[R, 5])+((Subscript[R, 1])+(Subscript[R, 3]))((Subscript[R, 2])+(Subscript[R, 4]))))//Simplify

In the case (Subscript[R, 5]) = 0 the result simplifies to

R =(((Subscript[R, 1]) (Subscript[R, 3]))/((Subscript[R, 1])+(Subscript[R, 3])))+(((Subscript[R, 2]) (Subscript[R, 4]))/((Subscript[R, 2])+(Subscript[R, 4])))

that can be obtained in a simple way. In the case (Subscript[R, 5]) = ∞ the result simplifies to

R =((((Subscript[R, 1])+(Subscript[R, 2]))((Subscript[R, 3])+(Subscript[R, 4])))/(((Subscript[R, 1])+(Subscript[R, 2]))+((Subscript[R, 3])+(Subscript[R, 4]))))

that also can be obtained elementarily.

Systems of linear equations in the matrix form

Systems of linear equations can be written in the matrix form

A.X==B,

where is the square matrix of the coefficients of the system of equations, X is the vector (List) of unknowns, and B is the vector (List) of right parts. In particular, for the problem of a rod on two supports described by the system of equations

(Subscript[F, 1])+(Subscript[F, 2])-M g==0    (* Sum of all forces is zero *)
(Subscript[F, 1])(Subscript[x, 1])+(Subscript[F, 2])(Subscript[x, 2])-M g L/2==0  (* Sum of all torques (here with respect to the left end) is zero *)

one has

A=(({{1, 1}, {(Subscript[x, 1]), (Subscript[x, 2])}}))
B={M g,M g L/2}   
X={(Subscript[F, 1]),(Subscript[F, 2])}

It looks natural to define

B=(({{M g}, {M g L/2}}))   
X=(({{(Subscript[F, 1])}, {(Subscript[F, 2])}}))

but, unfortunately, the current version of Mathematica understands this as Lists with two indices (Lists of Lists) that are inappropriate for us. Please, rerun now the previous definitions!

One can see that technically matrix is a List with two indices, whereas a vector is a list with two indices. The equations for the rod are

A.X==B

that is the same as equations written above.

The advantage of the matrix form is that properties of matrices are very well investigated and many operations can be performed on matrices as the whole rather than on their components.There are functions of matrices etc. The same pertains to vectors. If the problem is put into the matrix-vector form (that is, in the form of Lists) , one says that the problem is vectorized. Vectorized problems are solved faster because there are special highly efficient algorithms for Lists.

The solution of the rod problem with Mathematica can be obtained as

Solve[A.X==B,X]//Simplify

Also one can solve the matrix equation simply by left-multiplying the right and left parts by the inverse A matrix

Inverse[A]

and using

Inverse[A].A==(({{1, 0}, {0, 1}}))//Simplify

or, in other words

Inverse[A].A==IdentityMatrix[2]//Simplify

This yields the solution

X=Inverse[A].B//Simplify

Note that the command

(A^(-1))

does not yield the inverse matrix, it simply inverts each element. All operations on matrices have special names. Similarly, the exponential of a matrix is given by MatrixExp[A] and not by Exp[A] (the latter exponentiates each element).

Systems of nonlinear equations

Physically, systems of nonlinear equations arise in finding the equilibrium state of a system corresponding to the minimal potential energy and in other problems. Systems of nonlinear equations can be solved by Mathematica by the same methods as nonlinear equations. 

Consider, as an illustration, two particles in one dimension, described by the coordinates x and y that repel each other with the quadratic potential (Subscript[U, 12])[x,y]=-(((x-y))^2). For each particle there is a confinement potential (Subscript[U, 1])[x]=(x^4) (and the same for y) that prevents the particles from going too far from each other. As it is seen from the plot below, the minima of the energy correspond to {x=±1,y=∓1}.

xMax=1.5;
U[x_,y_]=-(((x-y))^2)+((x^4)+(y^4));
Plot3D[U[x,y],{x,-xMax,xMax},{y,-xMax,xMax}]

Extrema of a function of two variables satisfy

(∂_x)U[x,y]==0,       (∂_y)U[x,y]==0

that form a system of two equations for two unknowns

(∂_x)U[x,y]==0
(∂_y)U[x,y]==0

This system of equations can be solved analytically. Adding and subtracting the equations, one obtains the equivalent system of equations

(x^3)+(y^3)==0
(x^3)-(y^3)-(x-y)==0      ⇒      (x-y)((x^2)+xy+(y^2)-1)==0

From the first equation one obtains y=-x. Substituting this into the second equation, one obtains the equation for x

x((x^2)-1)==0

that solves to

x=0,±1.

Here x=0 is a local maximum while x=±1 are minima.

Solution of this problem with Mathematica that is good for polynomial functions

Solve[{(∂_x)U[x,y]==0,(∂_y)U[x,y]==0}]

yields all extrema, including unphysical complex solutions. Or, numerically,

NSolve[{(∂_x)U[x,y]==0,(∂_y)U[x,y]==0}]

FindRoot can find solve systems of transcedental equations, if the initial points are specified and are sufficiently close to the actual roots

FindRoot[{(∂_x)U[x,y]==0,(∂_y)U[x,y]==0},{{x,0.5},{y,-0.5}}]

Another example. Three beads electrically charged with positive charges (Subscript[Q, 1]), (Subscript[Q, 2]), and (Subscript[Q, 3]) can freely glide on a circular contour of radius R. Find their equilibrium positions. 

Solution. The potential energy of the system is a sum of Coulomb potential energies of point charges,

U[(Subscript[r, 1]),(Subscript[r, 2]),(Subscript[r, 3])] =(((Subscript[Q, 1])(Subscript[Q, 2]))/(|(Subscript[r, 1])-(Subscript[r, 2])|))+(((Subscript[Q, 2])(Subscript[Q, 3]))/(|(Subscript[r, 2])-(Subscript[r, 3])|))+(((Subscript[Q, 3])(Subscript[Q, 1]))/(|(Subscript[r, 3])-(Subscript[r, 1])|)).

This potential energy is invariant with respect to the rotation of the system of charges as the whole. Thus one can fix the position of (say) the third bead as {(Subscript[x, 3]),(Subscript[y, 3])}={R,0} (that is, (Subscript[ϕ, 3]) = 0) and describe the coordinates of the other two beads in the polar coordinate system

{(Subscript[x, 1]),(Subscript[y, 1])}={R Cos[(Subscript[ϕ, 1])],R Sin[(Subscript[ϕ, 1])]}
{(Subscript[x, 2]),(Subscript[y, 2])}={R Cos[(Subscript[ϕ, 2])],R Sin[(Subscript[ϕ, 2])]}
{(Subscript[x, 3]),(Subscript[y, 3])}={R,0}

The potential energy of the system can be rewritten as

U[(Subscript[ϕ, 1]),(Subscript[ϕ, 2])] =(((Subscript[Q, 1])(Subscript[Q, 2]))/(R(Sqrt[((((Subscript[x, 1])-(Subscript[x, 2])))^2)+((((Subscript[y, 1])-(Subscript[y, 2])))^2)])))+(((Subscript[Q, 2])(Subscript[Q, 3]))/(R(Sqrt[((((Subscript[x, 2])-(Subscript[x, 3])))^2)+((((Subscript[y, 2])-(Subscript[y, 3])))^2)])))+(((Subscript[Q, 3])(Subscript[Q, 1]))/(R(Sqrt[((((Subscript[x, 3])-(Subscript[x, 1])))^2)+((((Subscript[y, 3])-(Subscript[y, 1])))^2)])))

or, in terms of the two angles (Subscript[ϕ, 1])and (Subscript[ϕ, 2])

U[(Subscript[ϕ, 1]),(Subscript[ϕ, 2])] =(((Subscript[Q, 1])(Subscript[Q, 2]))/(R(Sqrt[(((Cos[(Subscript[ϕ, 1])]-Cos[(Subscript[ϕ, 2])]))^2)+(((Sin[(Subscript[ϕ, 1])]-Sin[(Subscript[ϕ, 2])]))^2)])))+(((Subscript[Q, 2])(Subscript[Q, 3]))/(R(Sqrt[(((Cos[(Subscript[ϕ, 2])]-1))^2)+(((Sin[(Subscript[ϕ, 2])]))^2)])))+(((Subscript[Q, 3])(Subscript[Q, 1]))/(R(Sqrt[(((1-Cos[(Subscript[ϕ, 1])]))^2)+(((Cos[(Subscript[ϕ, 1])]))^2)]))).

Now the potential energy can be plotted (Sorry, cannot define it as a function of (Subscript[ϕ, 1])and (Subscript[ϕ, 2]))

U[ϕ1_,ϕ2_] :=(((Subscript[Q, 1])(Subscript[Q, 2]))/(R(Sqrt[(((Cos[ϕ1]-Cos[ϕ2]))^2)+(((Sin[ϕ1]-Sin[ϕ2]))^2)])))+(((Subscript[Q, 2])(Subscript[Q, 3]))/(R(Sqrt[(((Cos[ϕ2]-1))^2)+(((Sin[ϕ2]))^2)])))+(((Subscript[Q, 3])(Subscript[Q, 1]))/(R(Sqrt[(((1-Cos[ϕ1]))^2)+(((Sin[ϕ1]))^2)])));
ParameterSet={R->1,(Subscript[Q, 1])->1,(Subscript[Q, 2])->1,(Subscript[Q, 3])->1};
Plot3D[U[ϕ1,ϕ2]/.ParameterSet, {ϕ1,0,2π},{ϕ2,0,2π},PlotRange->{0,3}]

Obviously, U diverges if two of three particles approach each other. For all charges being the same, as plotted, the minima correspond to {(Subscript[ϕ, 1])= 2π/3, (Subscript[ϕ, 2])= 4π/3} and {(Subscript[ϕ, 1])= 4π/3, (Subscript[ϕ, 2])= 2π/3}. It is interesting that Solve cannot find this obvious solution. Since there are no local maxima and saddle points, FindRoot works without problems, including the cases of different charges where it is impossible to find their equilibrium positions analytically. At the end we plot equilibrium positions of the beads.

ParameterSet={R->1,(Subscript[Q, 1])->0.3,(Subscript[Q, 2])->.1,(Subscript[Q, 3])->1};
MyRoots=FindRoot[{((∂_(ϕ1))U[ϕ1,ϕ2]/.ParameterSet)==0,((∂_(ϕ2))U[ϕ1,ϕ2]/.ParameterSet)==0},{{ϕ1,1},{ϕ2,2}}];

RR=R/.ParameterSet;
(Subscript[ϕ, 1])=ϕ1/.MyRoots;
(Subscript[ϕ, 2])=ϕ2/.MyRoots;
Graphics[{
Circle[{0,0},RR],
{Red,Disk[{RR Cos[(Subscript[ϕ, 1])],RR Sin[(Subscript[ϕ, 1])]},0.03]},
{Green,Disk[{RR Cos[(Subscript[ϕ, 2])],RR Sin[(Subscript[ϕ, 2])]},0.03]},
{Blue,Disk[{RR ,0},0.03]},
Disk[{0,0},0.01],
{Dashed,Line[{{RR ,0},{-RR ,0}}]},
{Dashed,Line[{{0 ,RR},{0,-RR}}]}
}]

The main algorithm used by FindRoot for systems of nonlinear equations is again the Newton method. Equations are linearized near the starting point and the resulting linear equations are solved to give the new starting point. If the starting point is close to the root, the convergence is very fast.