How Do I Solve Equations with Mathematica?

Polynomial Equations

Mathematica provides two main functions for solving polynomial equations. They are Solve and NSolve. Solve is used to algebraically solve an equation or set of equations. However, it is not always possible to solve an equation algebraically and there are times when a numerical solution is desired instead. In either case, NSolve can be used to find a numeric solution to an equation or set of equations.

Algebraic Solutions

Solve is the Mathematica function used for symbolically solving a polynomial equation or set of equations. Its syntax is Solve[eqns, vars], where eqns is your equation or set of equations and vars are the variable(s) in the equation(s). The equations are written in the form of lefthandside == righthandside. For example, the equation a (x^2)  +  b x  +  c  =  0 would be entered as

a (x^2)  +  b x  +  c  ==  0

Be sure to use the double equal sign when entering equations in Mathematica. For a more detailed explanation on using = versus ==, please see section 1.5.5 of the Mathematica book.

Single Equation

Let's solve the previously mentioned equation  a (x^2)  +  b x  +  c  =  0 for x.

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

Special Topic: Algebraic Numbers

The following sections give a high level introduction to the way Mathematica handles algebraic numbers. It is not essential to understanding the Mathematica equation solvers and is just included for interested readers. 
The text is adapted from a talk presented by Robert Knapp at the Japanese Society for Synmbolic and Algebraic Computing.

For higher order polynomials, it is not always possible to represent the exact solution in closed form. Currently, the highest polynomial that a closed form solution has been found for is the Quinoid or 5th order polynomial. However, the expression represents an exact mathematical quantity, which is based on the abstraction of algebraic numbers, represented in Mathematica as Root objects.

The method Mathematica uses internally to calculate roots of polynomials is the well established Jenkins-Traub method.  It works by finding successive roots with what is essentially a numerically stable variant of the Newton-Raphson method. Then the polynomial is deflated by dividing out by a factor involving each root that is found.  This is the engine behind algebraic numbers.  What makes it possible to have an implicit representation which is reasonably efficient are mathematical theorems which allow roots to be isolated and Mathematica's significance arithmetic.

Lets take a closer look at the roots of a higher order equation.

Solve[(x^5) - (x^3) +  σ x  +  1  ==  0, x]

gives the solution to the equation.  This is given in terms of replacement rules so that you can see how the variables relate to the solution.  In this case it indicates that there are five possible solutions.  In cases where you just want to work with the roots, you can extract them into a list.

polyroots  =  x  /.  %

A plot of the root positions in the complex plane as a function of the parameter value helps to visualize the numbers implicitly represented here. First, an auxiliary package is loaded in order to highlight the roots with different colors.

Needs["Graphics`Colors`"];

ParametricPlot[
        Evaluate[{Re[#], Im[#]} &  /@  roots], {σ, -2, 2}, 
        PlotStyle -> {Red,  Blue,  Green,  Black,  Purple}];

It is interesting to note that  two roots come together on the real axis, at the point σ = -1.  The lines going across the plot are an artifact of the canonical order Mathematica gives to the roots.  You may be asking: how does Mathematica distinguish between and order the roots?  The process is an elegant example of numerical and symbolic mathematics working together.  We will work with the last root to illustrate this and to make things a bit more convenient we define the root as a function of the parameter:

lastroot[σ_]  =  Last[roots];

When the parameter σ is given a numerical value, Mathematica does the computations necessary to determine an approximate root which is precise enough to be isolated from the other roots of the polynomial.

lastroot[0]

First, roots approximations at machine precision are found using the Jenkins-Traub algorithm.  The isolation is possible because for a given root approximation it is possible to find an enclosing circle in which an exact root of the polynomial is guaranteed to lie.  Isolation amounts to making sure that none of the circles around the approximate values overlap. Optionally Mathematica allows you to use an exact root isolation algorithm which is more complicated and time consuming.  See [Strzebonski, 1996] for details on both methods of isolation.  When σ = 0, the roots are reasonably far apart from each other, so that 16 digits is sufficient to ensure isolation.  Once the roots are isolated, they are ordered based on a canonical ordering.  Real roots are ordered first along the real line, then complex ones based on complex absolute value.  This is why you see the lines across the plot above at σ = -1.  The pairs of roots switch to maintain the canonical ordering.

When σ = -1, the polynomial is reducible in terms of radicals, and Mathematica recognizes this, so

lastroot[-1]

is given in its simplest from.  However, for σ sightly greater than -1, the polynomial cannot be reduced, and so

lastroot[-1 + (10^(-30))]

is given in terms of a root object. Notice that the polynomial is recalculated to be the minimal polynomial with integer coefficients.  Internally Mathematica recognizes that the two roots are very close together and adapts to using higher precision (in this case 32 digits) to do the root isolation for these.  Then, when you request a numerical value,

N[%]

the positive imaginary part is indicated correctly.

Hopefully this example has given you a rough idea of how the fusion of implicit representations with numerical approximations using adaptive precision, yields a powerful computational tool.  Mathematica provides a framework where objects like algebraic numbers fit naturally.  You can similarly use Mathematica to set up implicit representations of your own objects by defining numeric and symbolic rules which apply to them.  Implicit representations can be thought of in some sense as a transparent interface to a combination of powerful numerical algorithms and mathematical information without overburdening a user with the implementation details.

Special Topic: Finding Special Solutions

Solve is set up to give you generic solutions to equations. It discards any solutions that exist only when special constraints between parameters are satisfied. However, the function Reduce will keep all the possible solutions to a set of equations, including those that require special conditions on parameters.

Reduce[a x^2 - b == 0, x]

In this case the output means the possible solutions are

(1) a = 0 and b=0
or
(2) x = -((Sqrt[b])/(Sqrt[a])) and a is not equal to 0
or
(3) x = ((Sqrt[b])/(Sqrt[a]))and a is not equal to 0

Simultaneous Equations

Solve is also used in algebraically solving sets of simultaneous equations. Let's solve the two equations (x^2)  +  a x + (cy^2) = 0 and y + x = 5 for x and y.

Solve[{(x^2) + a x + c (y^2) == 0,  y + x ==  5}, {x, 
        y}]

Note that in the syntax above, braces surround both equations and both variables. Another way of entering the equations would be to assign the set of equations a name. For example,

\(\(eqns = {x\^2 + y\^2 + z\^2 == 20, 
	\ x + z == \ 
            2, 
	\ y + x == 0};\)\)

Solve[eqns,  {x, y, z}]

More on Solve.

Numeric Solutions

There is more than one way to get a numeric solution to an equation. For example, to numerically solve the equation (x^2) + 2 x  - 7 = 0, we could do either of the following.

First, we solve the equation algebraically and then use N to get a numeric answer.

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

N[ % ]

Or, we could use NSolve to numerically solve the equation from the start.

NSolve[(x^2) + 2 x - 7 == 0, x]

Here's an example of solving a more complicated polynomial equation using NSolve.

poly = 
        3 + 3  x - 7  (x^2) - (x^3) + 2  (x^4) + 3  (x^7) - 3  (x^8) - (x^9) + 
          (x^10);

NSolve[poly == 0, x]

NSolve can also be used to numerically solve systems of equations. The syntax is the same as for Solve.

NSolve[{x + y == 2, x - 3 y + z == 3, x - y + z == 0},
         {x, y, z}]

More on NSolve.

Special Topic: Large Systems of Polynomial Equations

In many engineering, financial and scientific calculations, researchers have to deal with large systems containing hundreds or thousands of equations. In those cases, linear algebra methods are usually more suited than Solve or NSolve. These methods are covered in Section 3.7.8 through Section 3.7.10 of The Mathematica Book. If you are dealing with sparse systems, please read the documentation for Developer`SparseLinearSolve.

Other Equations

In addition to being able to solve purely algebraic equations, Mathematica can also solve some equations involving other functions. Solve can be used in solving radical equations, equations involving trigonometric or hyperbolic functions and their inverses as well as equations involving exponentials and logarithms. As a reminder its syntax is NSolve[eqns, vars], where eqns is your equation or set of equations and vars are the variable(s) in the equation(s).

Algebraic Solutions

Radical Equations

Solve[(1/x) + (1/(Sqrt[x])) == 1, x]

Note that in radical equations, Solve discards parasite solutions. To see all candidate solutions, including parasites, set VerifySolutions to False.

Solve[(1/x) + (1/(Sqrt[x])) == 1, VerifySolutions -> False]

Equations involving trigonometric or hyperbolic functions and their inverses

Solve[{Cos[x + α + (1/5)] == 5, 
        Cos[-x + α + (2/5)] == 6}, {α, x}]

Equations involving exponentials and logarithms

Solve[(((1 + (x^3)))^(a + x)) == (E^(a + x)), 
      x]

More on Solve.

Numeric Solutions

If your equations involve only linear functions or polynomials, then you can use NSolve to get numerical approximations to all the solutions. However, when your equations involve more complicated functions, there is in general no systematic procedure for finding all solutions, even numerically. In such cases, you can use FindRoot to search for solutions. The basic syntax for FindRoot is FindRoot[eqn, {x, (Subscript[x, 0])}] where eqn is the equation you are solving and (Subscript[x, 0]) is the value around which FindRoot starts its search.

FindRoot[ 3 Cos[x] == Log[x], {x, 1} ]

More on FindRoot.

Differential Equations

Please see the How do I solve differential equations in Mathematica? how to.

More Information

The Mathematica book and the Built-in Functions section of the Help Browser contain much more information on solving equations in Mathematica. The following are some suggested sections to take a look at if you want to know more about solving equations.

The Mathematica Book Sections
1.5.5
1.5.7
1.5.8
1.6.3
1.6.4
3.4
3.5.10
3.9.5
3.9.6
3.9.7
3.9.8

Built-in Functions in Help Browser
Solve
NSolve
DSolve
NDSolve
FindRoot