Tactics for solving equations in Mathematica

Daniel Lichtblau
Wolfram Research, Inc.
danl@wolfram.com
Original version: June 1998
Most recent update: January 2005

Abstract

We will discuss equation solving via Solve and other commands. We will outline the workings of Solve, give simple examples of its use, and provide pointers to related functions and when/how to use them. We will also show ways in which to tackle problems that are not immediately tractable for Solve.

Initialization

Off[General::spell]; Off[General::spell1];

Functions in Mathematica for solving equations and systems of equations

The general purpose algebraic equation solver engine is Solve. For linear systems one might instead use matrix notation and call LinearSolve or  RowReduce (an example of the difference will be shown below). For large sparse systems with machine number real or constant coefficients, there is Developer`SparseLinearSolve. For solving linear systems repeatedly with the same matrix of coefficients one might use LUDecomposition followed (repeatedly) by LUBackSubstitution. For approximate solutions to algebraic systems that have finitely many solutions there is NSolve. For approximate solutions near a given point one can use FindRoot or, in some cases, FindMinimum. It is not feasible to discuss all these in detail. I will show examples of several and for the rest we will make do with general remarks.

How Solve works

I will begin with an abbreviated description of what Solve does. First we look for special cases such as one polynomial equation in one unknown, or linear equations in the given variables. Then we try preprocessing to transform various sorts of transcendentals such as logs or exponentials. We also will attempt factorization and other methods of regrouping to try to obtain simpler equations. After this we put the equations into disjunctive-normal-form (Or's of And's). After this we make the equations into polynomials in terms of internal variables. Transcendentals, radicals, and denominators will all become new variables, with appropriate auxiliary polynomials that describe relations between them and other variables. For example, Sqrt[x-1] gives rise to an equation newvar^2-(x-1). Cos[x+y] will get a new variable that in essence "inverts" it (using InverseFunction) so that we can get at it's innards.

At this point we sort variables and form a lexicographic Gröbner basis for the polynomials with respect to the ordered variables. We then do further work to try to determine if the system has finitely many solutions (did we solve for all variables?) because there are some simplifications in this case. After this we perform root extraction. This itself is a strenuous procedure. It may try factoring, decomposition, Cardano formulas (with attempts to simplify radicals used therein), and algebraic number (Root[...]) technology. Roots obtained for the last variable will be back-substituted into polynomials in the basis in the remaining variables, we get roots for the next-to-last variable, back substitute, etc.

Once we have found solutions we need to check for circular dependencies (that is, did we solve for a variable in terms of itself?) in cases where we had transcendentals in the input. Under some circumstances we may also need to attempt numerical verification, in case non-polynomial functions were in the input (say radicals or transcendentals) because in such cases we may get parasite solutions from the process of forming polynomials of the input.

Types of equations Solve can frequently handle

Below are several classes of problem that can frequently be handled simply by calling Solve.

(1) Polynomial systems (it is designed for this).

(2) More generally, algebraic systems (that is, we allow radicals). Here is an example.

Solve[(a^2) - (Sqrt[a]) + (a^(2/3)) == 5, a]

(3) Systems that are polynomials in trigs, exponentials, or (sometimes) logs. Also, equations that mix a power and an exponential or logarithm can frequently be solved in terms of the ProductLog function. Here are a few examples, taken from benchmark tests. The last example below is a power-and-logarithm problem in disguise (take logs of both sides). Further down I show a related example, where we use an infinite ladder of exponentials in the variable of interest.

Solve[{(E^(x + 3 y - 2)) == 7, 
        (3^(2 x - y + 4)) == 2}, {x, y}]

Solve[
      -Log[3] + (1/2) Log[2 + (y/3)] - (1/2) Log[(y/3)] == 
        (1/2) (-I) π, y]

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

Solve[(a^(2 x)) - (a^(x/2)) + (a^(x/3)) == 2, x]

Solve[((((x^2) - 5 x + 5))^((x^2) - 9 x + 20)) == 1, x]

Solve[(x^x) == E, x]

(4) Systems that can be transformed into amenable problems of the above type by adroit use of Collect, Together, Factor, etc. A fair amount of transformative (is that a word?) preprocessing is built into Solve; we try to look for excessive leaf count growth in an attempt to avoid hanging due to such transformations.

Solve is by no means perfect in regards to (3) or (4), but it is tolerably good. Several examples may be found in the "further examples" for Solve in the on-line documentation.

(5) Simple equations involving Abs, Min, etc. Anything too hard for Solve might be attempted with Algebra`InequalitySolve`. Actually this is often preferable for the simple reason that it can represent solution sets in regions of real space, unlike Solve.  Examples:

Solve[Abs[x - 1]  ==  Min[x + 2,  4],  x]

Solve[Abs[x - 1] + Abs[x + 2] == 3, x] (* gives garbage result *)

<< Algebra`InequalitySolve`InequalitySolve[Abs[x - 1] + Abs[x + 2] == 3, x]

InequalitySolve[Abs[x - 1]  ==  Min[x + 2,  4],  x]

We do not pursue inequalities further here. This is the topic of another talk and Mathematica notebook.

Tactics for solving difficult systems of equations

(1) Subtract right-hand-sides from left-hand-sides, apply Numerator[Together[#]]&, set result equal to zero. This will ignore difficulties at vanishing denominators but is often more efficient than built-in methods. The built-in code will create extra variables and equations for every denominator, in effect as newvar*denom==1 (this enforces that denom cannot vanish). These extra equations can drastically slow the solve process.

(2) For strenuous nonlinear polynomial systems, precompute GroebnerBasis[polynomials,vars], set the result to zero, and call Solve on that. This is often effective because Solve is not yet overhauled to use the better Gröbner code. More importantly, you can play with the variable ordering and might find one that is better than the ordering Solve chooses according to its heuristics.

(3) Use FindRoot. Many problems do not have "nice" algebraic solutions, and FindRoot will use Newton's method and relatives (e.g. secant method) to find an approximate solution, given a reasonable starting point. Caveat: FindRoot requires same number of equations as unknowns. If for some reason the system is underdetermined then you may want to "specialize" some variables, that is, choose particular values for them, and then solve for the rest. Alternatively you might try FindMinimum. Below is a simple example. We want a solution to one equation, x^2+y^2==1, in two variables x and y. If we only require a real solution we can rephrase it as a FindMinimum problem by squaring the (left hand side) - (right hand side). If we had more equations we would sum all such squares.

FindMinimum[(x^2 + y^2 - 1)^2, {x, 3}, {y, -4}]

(4) If you have algebraic equations (nothing transcendental), only finitely many solutions, want all of them, and will settle for approximations, try NSolve. Our new NSolve code is generally alot faster than Solve for such systems.

(5) If you have a symbolic linear system you can expect intermediate coefficient swell as the culprit of hanging or huge solution sets. For such systems Solve will directly call RowReduce (not LinearSolve, more on this below). Hence tactics to improve matters might involve setting the Method option of RowReduce to either CofactorExpansion or OneStepRowReduction. The default of Automatic will generally select DivisionFreeRowReduction because it is often fastest on "small" examples and because I know of no heuristics to decide between methods for symbolic linear algebra. For systems that have symbolic entries that are low-degree polynomial or rational functions in one variable you might go so far as to code something to interpolate the solutions. The result is often remarkably smaller and more rapidly obtained than what you might get using built-in methods of row reduction. See The Mathematica Journal, Vol. 7 #1, pp 26-27 for an example of how one might do this.

(6) Sometimes conversion from trigs to exponentials, or expansion of trigs, can be of help. For example, the following was posed in the news group comp.soft-sys.math.mathematica.

eqns  =  
      {-L + C*x + (2*R)/Tan[y/2], 
        -x + 2*((V - R - R/Cos[y])*Tan[y])};

Timing[sol1 = Solve[Together[TrigToExp/@eqns] == 0,  {x, y}]; ]

An alternative that works better is to expand the trigs, first changing y to 2*z to get rid of the denominator inside one of the trigs. We get a result with small LeafCount quite fast.

eqs = 
        TrigExpand/@
          {-L + C*x + (2*R)/Tan[z], 
            -x + 2*((V - R - R/Cos[2*z])*Tan[2*z])};Timing[
      sol = Solve[eqs == 0, {x, z}] /. 
          (z -> a_) :> y -> 2*a; ]

\(us\  = \ 
      \((u*\((1\  + \ 2*v\  + \ h*t*\((\(-1\)\  + \ θ)\))\)\  + 
	
            v*\((1\  + \ h*t*\((\(-1\)\  + \ θ)\)\  - \ t*θ)\))
          \)/
	\((
          2\  + \ u\  + \ v\  + \ 
            t*\((\(-2\)\  + \ 
                  h*\((u\  + \ v)\)*\((\(-1\)\  + \ θ)\)\  - \ 
                  v*θ)\))\)\)

Working through some difficult systems by hand

We will discuss a few systems of equations that were posted to the Usenet news group comp.soft-sys.math.mathematica. These use methods discussed above for making a system more tractable, analyzing it, and the like. Much of the methodology shown below comes from responses posted by the author to that news group.

example (i)

We start with the following set of equations.

exprs = {x1 x2 - x3 x4 - 100, (3/x2) - (x1^2) - (1/(10^3)), 
        (x4/x1) + x2 x3 - (2/(10^6)), x1 x4 - (x3/x2) - (5/(10^13))};

The poster was using FindRoot to look for solutions, and  noted that machine precision  did not seem to give adequate results. Two approaches worth trying are to use Solve or NSolve directly, or to use FindRoot with some sensible starting values. It happens that the system above has no real solutions, hence giving real starting values is hopeless because FindRoot uses variations of Newton's method and will never leave the real axes in its search for a result.

First we show that Solve can do this directly.

\(vars\  = \ {x1, x2, x3, x4}; 
\ \ sol\  = \ Solve[exprs == 0, \ vars]; 
    
	nsol1\  = \ N[sol, \ 100]; 
\ \ \ exprs\  /. \ nsol1\  // \ Chop\)

This also shows that sometimes the default methods are not quite adequate. If we do

SetOptions[Roots,  ExactRootIsolation -> True]

and then redo the Solve and N[...] all will be better when we plug the solutions into the original expressions.

We now use FindRoot to obtain a solution. The method below will work for some random starting values, but not for others. All the same, it gives a way to get solutions without undue labor.

\(\ max\  = \ 50; 
\ \ \ 
    start\  := \ Table[Random[Complex, max + I*max, 100], \ {4}]\)\(\ \ \ 
    FindRoot[Evaluate[Thread[exprs == 0]], 
\ \ \ \ \ \ \ \ \ \ \ 
      Evaluate[Sequence\ @@\ Thread[{vars, start}]], 
\ \ \ \ \ \ \ \ \ \ \ 
      WorkingPrecision -> 100, \ AccuracyGoal -> 20, 
\ \ \ \ \ \ \ \ \ \ \ 
      MaxIterations -> 300]\)

One comment is that, given the spread in scales of the input coefficients, machine precision will most likely be inadequate for getting results. Hence the desire to meddle with the WorkingPrecision, AccuracyGoal, and MaxIterations options of FindRoot.

example (ii)

The following system is a bit tricky. The user is interested in positive solutions, should any exist. It looks like the FindRoot method from example (i) above might work for it, but this proves not to be the case.

polys = {(369664/981) - (1/(x2 (x3 + x4))), 
        -(6513786578938559/9223372036854775808) + x1 x4 + x1 x5 + 
          x6 x7, -(3858337547701687/9444732965739290427392) + 
          x1 x2 x4 x5 + x1 x4 x6 x7 + x1 x5 x6 x7, 
        -(1/97280000000) + x1 x2 x4 x5 x6 x7, 
        (981/369664) - x2 x3 - x2 x4, 
        (77/924160000) - x1 x2 x3 x4 - x1 x2 x3 x5 - x2 x3 x4 x6 - 
          x2 x3 x6 x7 - x2 x4 x6 x7, 
        (1/1848320000000) - x1 x2 x3 x4 x6 x7 - x1 x2 x3 x5 x6 x7}; 
    vars = {x1, x2, x3, x4, x5, x6, x7};

FindRoot gives failure-to-converge warning and a bad residual. When the
 MaxIterations option is upped considerably we still have this problem. We will see why presently. First we prepare to use exact solving methods, by clearing denominators, extracting numerators, and forming a Gröbner basis.

\(p2\  = \ Numerator[Together[polys]]; 

    gb\  = \ GroebnerBasis[p2, \ vars]; 
Length[gb]\)

This reveals the source of the difficulty in using FindRoot. We have an underdetermined system. The solution set is one dimensional, and Newton's method fails because the Jacobian is singular. All the same, for this particular system we can do alot more.

gb[[1]]

We see that the first polynomial can be made univariate by substituting z for x6*x7.

pol1 = 
      Last[PolynomialReduce[gb[[1]], 
          x6 x7 - z, {x6, x7, z}]];

We now find the roots of this polynomial.

rz  =  NRoots[pol1 == 0,  z,  50]

We have a positive root, which means we might have solutions for the original system for which x6 and x7 are positive, and their product must then be approximately 0.0000263158. We can now play with FindRoot a bit more, giving these new polynomials and adjoining a specific value for x7, say, to have #equations = #unknowns.

I chose to set x7 to one. By trial and error I found I had bad luck with MaxIterations set to 100, but better results at 500. Caveat: this seems to depend heavily on having "lucky" random starting values.

\(max\  = \ 50; 

    start\  := \ Table[Random[Real, max, 100], \ {Length[vars]}]\)\(\(rt\  = \ 
      FindRoot[Evaluate[Thread[Prepend[gb, x7 - 1] == 0]], 

        \ \ \ \ \ \ \ \ \ \ \ Evaluate[Sequence\ @@\ Thread[{vars, start}]], 
        
\ \ \ \ \ \ \ \ \ \ \ WorkingPrecision -> 100, \ 
        AccuracyGoal -> 20, 
\ \ \ \ \ \ \ \ \ \ \ MaxIterations -> 400]; 
    \)\)

polys  /.  rt

We got real roots, and small residuals in the bargain. But the roots are not all positive. By looking hard at gb, and keeping in mind what the product x6*x7 must be, we can prove that there is no nonnegative solution to the system.

Last/@
      PolynomialReduce[gb, 
        x6 x7 - rz[[1, 2]]]

As expected, we get a small value for the first polynomial, because it is supposed to vanish for that value of x6*x7. The rest give polynomials that are each linear in x7 and in one other variable, in effect solving for that variable in terms of x7. We see, moreover, that all coefficients are positive. So after dividing each polymonial by its x7 coefficient we have eqns of the form positive*xj+x7==0. This gives a negative solution for each xj for 1≤j≤5. Conclusion as stated above: there is no root of the original system with positive-only values. Which agrees with my experience using FindRoot on gb several times in the manner shown above.

example (iii)

Below is a difficult problem with many parameters that was posed in the news group. We show explicitly the various tactics, discussed above, that were attempted by the user and others to make it tractable for solving.

us = 
      ((u (1 + 2 v + h t (-1 + θ)) + 
          v (1 + h t (-1 + θ) - t θ))/(2 + 
          u + v + t 
            (-2 + h (u + v) (-1 + θ) - v θ)
              )); vs = 
      ((u (1 + 2 v + k s (-1 + θ)) + 
          v (1 + k s (-1 + θ) - s θ))/(2 + 
          u + v + s 
            (-2 + k (u + v) (-1 + θ) - v θ)
              ));

The variables of interest are u and v. The rest are constants such that
(0 ≤ k ≤ 1, 0 ≤ h ≤ 1, -1 < s < 1, -1 < t < 1). 
If need be we can just handle the case of most interest, where h==k.

The poster first tried one of the recommendations above, that is, to clear denominators.

\(equations\  = {us == u, vs == v}; 

    polys = Map[#[\([1]\)] - #[\([2]\)]&, equations]; 

    newpolys\  = \ Map[Numerator[Together[#]]&, polys]\)

As it happens, this was not sufficient to yield a system tractable to Solve. We will explicitly compute a Gröbner basis, setting the CoefficientDomain option so that we treat the parameter variables as coefficients rahter than as variables.

\(vars\  = \ {u, v}; 

    gb\  = \ GroebnerBasis[newpolys, \ vars, 
	
        CoefficientDomain -> RationalFunctions]; \)

We now check the result to see if it is usable.

Length[gb]LeafCount[gb]Map[Exponent[#, u]&,  gb]Map[Exponent[#, v]&,  gb]

We have a univariate quartic polynomial in v, and a linear polynomial in u with coefficients that are at worst cubics in v. This means in principle we can solve for v using radicals or Root objects, then back substitute to solve for u. Now with polynomials of this size you do not want to attempt radical solutions. But even if we inhibit Roots from using the complicated radical formulas, it still appears to hang. So what can we do? We might substitute values for the parameters to reduce the complexity of the coefficients. First we will simplify the polynomials. The best result I obtained was as below.

\(collectedgb\  = \ Map[Collect[#, vars, Simplify]&, \ gb]; 

    LeafCount[collectedgb]\)

This may still be too complicated to handle in Solve or Roots. It is now at least in a form where one can look at the polynomials, look for reasonable substitutions, parameter specializations, and the like.

Reasons and tactics for numeric solving

Often Solve either cannot handle a problem due to transcendental dependencies, or else there are numeric methods that might be preferable. In the previous section the exact solution tactics were emphasized. Here we work primarily with numeric methods e.g. FindRoot and NSolve.

example (i)

For example, a question that came to the Usenet news group comp.soft-sys.math.mathematica was as follows. "I need to plot the root of a transcendental equation Sin[x]==x/a, for a from 2 to 3. How might this be done?" A simple answer uses FindRoot.

f[a_]  :=  FindRoot[Sin[x] == x/a,  {x, 2}]   Plot[x  /. f[a],  {a, 2, 3}];

example (ii)

Another reason to prefer a numeric solver might be when one particular solution is desired. Below is an example, also from comp.soft-sys.math.mathematica. The user wants only positive solutions for a given that c is an integer around, say, 100.

P[c_, a_] := 
      ((a^c)/(c! (∑_(i = 0)^c) ((a^i)/(i!))))FindRoot[P[100, a] == (1/50), {a, 10}, WorkingPrecision -> 50]

example (iii)

A related question, also posed to the Mathematica users' news group, was as follows. I have a  family of functions vv[m]:[0,1]-->R, parametrized by m, which may take any real value in the interval from one to infinity.  For each m, vv[m] is defined  to be(1-b-Exp[-b*m])/b for 0<b<=1, and m-1 when b is zero. I want to define a parametrized inverse bb[m] such that bb[m][vv[m][b]]==b and vv[m][bb[m][v]]==v.  Note that for each m, the function vv[m] is strictly monotonic in b, hence the desired functions bb[m] are all well-defined. How might I obtain this parametrized inverse in Mathematica?"

One can do this directly in Mathematica, but as we show below the result may not be quite what is wanted.

eq = vv == ((1 - b - (E^(m b)))/b); 
    inv = Solve[eq, b]

The ProductLog is the multi-valued function that inverts y==x*Exp[x] to obtain x as a function of y (it is sometimes called Lambert's W function in the literature). We actually have a family of solutions of the form
{{b -> (m/(1 + vv) - ProductLog[k, (E^(m/(1 + vv))*m)/(1 + vv)])/m}}
indexed by integers k. Graphing the ProductLog function for the zero (default) and minus one branches shows real values near the origin. One can see this as below.

\(Off[Plot::plnr]; 
Plot[ProductLog[x], \ {x, \(-1\), 1}]\ ; 

    Plot[ProductLog[\(-1\), x], \ {x, \(-1\), 1}]; \)

Warning messages (which I supressed) and the plots themselves indicate where each branch runs into trouble.

We now redo this numerically, and later we will compare results.

invert[m_] := 
      FindRoot[#1 == ((1 - x - (E^(m x)))/x), {x,  .7}]&v[m_] := ((1 - b - (E^(m b)))/b)

Now let's try to see what might be a reasonable value to invert.

Plot[v[3. ],  {b,  .1,  1}];

A look at the plot reveals that in the relevant range 0<b<=1 we take on
values in the range around -4 to -20. We will invert at -10.

invert[3. ][-10]

Now let us see if this is consistant with the result from Solve.

inv  /.  {m -> 3. ,  vv -> -10. }

Not even close. We will change the Solve result by hand, using the minus one (rather than default, zero) branch of ProductLog, and see if this helps.

inv2 = {
        b -> 
          (((m/(1 + vv)) - 
              ProductLog[-1, 
                (((E^((m/(1 + vv)))) m)/(1 + vv))])/m)};

inv2  /.   {m -> 3. ,  vv -> -10. }

So the branch used by Solve was not the one we really wanted after all. Really it was just luck (good or bad is for you to determine) that Solve could give a result at all. The root-finding method seems to be the better way to go for this sort of problem.

example (iv)

Frequently one has a system of polynomial equations with finitely many solutions. In such cases it is often most efficient to use NSolve rather than Solve. Below are some examples from our tests.

polys = {15 (b^2) c (d^2) + 6 (b^2) (c^3) + 21 (b^2) (c^2) d - 
          144 b c - 8 b (c^2) e - 28 b c d e - 648 b d + 
          b (d^2) e + 9 (b^2) (d^3) - 120, 
        30 (b^2) (c^3) d - 32 c d (e^2) - 720 b c d - 24 b (c^3) e - 
          432 b (c^2) + 576 c e - 576 d e + 16 b c (d^2) e + 
          16 (d^2) (e^2) + 16 (c^2) (e^2) + 9 (b^2) (c^4) + 5184 + 
          39 (b^2) (c^2) (d^2) + 18 (b^2) c (d^3) - 432 b (d^2) + 
          24 b (d^3) e - 16 b (c^2) d e - 240 c, 
        216 b c d - 162 b (d^2) - 81 b (c^2) + 5184 + 1008 c e - 
          1008 d e + 15 b (c^2) d e - 15 b (c^3) e - 80 c d (e^2) + 
          40 (d^2) (e^2) + 40 (c^2) (e^2), 
        261 + 4 b c d - 3 b (d^2) - 4 b (c^2) + 22 c e - 
          22 d e}; 
vars = {b, c, d, e}; 
Timing[
      solns200 = 
        NSolve[polys == 0, vars, WorkingPrecision -> 200]; ]

Apply[Plus, Abs[polys /. solns200]]

c1 = {x - (Sqrt[(t^2) - 1]), y - ((t^3) + t - 4)}; 
c2 = {x - ((s^2) + s + 5), y - (Sqrt[(s^2) + 7 s - 2])}; 
NSolve[Join[c1, c2], {x, y}, {s, t}]

Comparison of methods for solving linear systems

Much of this section may be found in the further examples of the on-line documentation. All the same perhaps it bears repeating in a basic notebook about equation solving in Mathematica.

Here is a simple LinearSolve example.

LinearSolve[{{1, 2}, {1, 3}}, {5, 8}]

You can also use Solve on the system.

Solve[{x + 2 y == 5, x + 3 y == 8}, {x, y}]

You get a similar result from RowReduce by creating an augmented matrix. Interpret the result by noting that the first column corresponds to the 
variable x and the second to the variable y; the solution is in the last 
column.

RowReduce[{{1, 2, 5}, {1, 3, 8}}]

For underdetermined systems Solve turns the excess variables into "
parameters" and gives general solutions in terms of those parameters.

Solve[{x + 3 y + z == 1, 2 x + y - 5 z == -2}, 
      {x, y, z}]

LinearSolve, on the other hand, gives a particular solution.

LinearSolve[{{1, 3, 1}, {2, 1, -5}}, {1, -2}]

RowReduce gives a similar result to Solve.

RowReduce[{{1, 3, 1, 1}, {2, 1, -5, -2}}]

This is interpreted as {x-16/5*z==-7/5,y+7/5*z==4/5} which indeed agrees with the solution from Solve. About the only way it might not agree is if we did not specify our variables to Solve and it ordered them in such a way that a different one was selected as the parameter.

This is a good place to mention LUDecomposition. For systems that are exactly determined (that is, n equations in n unknowns with exactly one solution), this can be a useful tool. Here is an example. We use an asymmetric variation of a 15×15 Hilbert matrix, with approximate entries.

h = N[Table[1/(2*i + j - 1), {i, 15}, {j, 15}], 30];

For our right-hand-side we will use a unit vector along the first coordinate.

e1 = Table[If[j == 1, 1, 0], {j, 15}];

We solve the system h.x==e1, first via LU techniques and then with 
LinearSolve.

lud = LUDecomposition[h];

As an aside I note that the third element returned is an approximation of the matrix condition number. Its logarithm (that is, scale) gives an indication of how many digits of precision one might expect to lose in solving linear equations with that matrix on the left-hand-side. Since machine arithmetic usually has about 16 places, something larger than, say, (10^15) indicates strong likelihood that machine arithmetic will return completely bogus results. The condition number of h is around (10^22), hence we must use bignum arithmetic.

Last[lud]

Back to solving.

sol1  =  LUBackSubstitution[lud, e1]

sol2  =  LinearSolve[h, e1]

sol1  -  sol2

Note that the solutions agree to at least 15 significant places in all coordinates.

Series solutions

Sometimes symbolic solving is not possible. In such cases a series solution might be attempted. In the example below, also from a benchmark test suite in the literature, we use a series "ansatz" to obtain a full solution. We wish to find constants {a,b,c,d} such that
a*Exp[b*x] - c*Exp[d*x] is equal to Sin[x]. We expand the difference as a power series to get polynomial equations, solve them, and this gives the correct result.

sols = 
      Solve[Series[
            Sin[x] - a (E^(b x)) - 
              c (E^(d x)), {x, 0, 10}] == 0, {a, b, c, d}]

Now check the result.

MapAll[ExpandAll[TrigToExp[#]] &, 
      Sin[x] == a Exp[b x] + c Exp[d x] /. sols]

Obtaining "full" solution sets to certain transcendental equations

We are asked from time to time whether Solve can find ALL solutions to e.g. Cos[x]==1. That is, users sometimes want a convenient way to denote solution sets that involve multi-valued inverses. Mathematica only gives principle inverses to e.g. trigonometric functions, so some work needs to be done to get at other potential solutions.

Something along these lines can often be done in the manner indicated below. For each solution we will generate a unique string prefixed by the character K. It is meant to stand for "any integer."

When we do this it pays to remove the Solve::ifun message. Moreover there are several caveats that are demonstrated in later examples.

arctrigs = {{ArcSin, 2}, {ArcCos, 2}, {ArcCsc, 2}, {ArcSec, 
            2}, {ArcTan, 1}, {ArcCot, 1}, {ArcSinh, 2*I}, {ArcCosh, 
            2*I}, {ArcCsch, 2*I}, {ArcSech, 2*I}, {ArcTanh, I}, {ArcCoth, 
            I}};
Do[arcTrigQ[arctrigs[[j, 1]]] = True;  
        mult[arctrigs[[j, 1]]] = arctrigs[[j, 2]], {j, 
          Length[arctrigs]}];
generalize[(f_)[x_]] /; TrueQ[arcTrigQ[f]] := 
      f[x] + mult[f]*Pi*Unique[K]
generalize[Log[x_]] := 
      Log[x] + 2*Pi*I*Unique[K]
generalize[ProductLog[x_]] := 
      ProductLog[Unique[K], x]
generalize[x___] := x
Off[Solve::ifun]

Here are some examples.

soln = 
      generalize //@ 
        Solve[(E^(2 x)) - 2 (E^x) + 2 == 0, 
          x]

Simplify[(E^(2 x)) - 2 (E^x) + 2 == 0 /. 
        soln, Element[Variables[x /. soln], Integers]]

eqns = {tt + aa == 46, tt + aa (E^(10 x)) == 39, 
        tt + aa (E^(20 x)) == 33}; 
    First[generalize//@Solve[eqns, {x, tt, aa}]]

Be warned that we still may not get all solutions in this manner. We will not get any if, for example, the non-generalized solution set is empty.

generalize//@
      Solve[{(E^(I x)) == I, 
          (y/(x - (π/2))) == 3}, {x, y}]

Moreover we fail to obtain generalized solutions in cases where an inverse trig evaluates to something numeric.

generalize //@ Solve[Cos[x] == 1,  x]

This is in contrast to the case where the inverse trig remains unevaluated.

generalize//@Solve[Cos[x] == 2,  x]

A simple way to obtain solutions for the first case is to keep it symbolic and substitute a numeric value after-the-fact. In this simple example we get two equivalent families of solutions.

(generalize //@ Solve[Cos[x] == a,  x])  /.  
      a -> 1

We can also get spurious solutions. In the example below, the "solution" for x is bad when the integer parameter K$xx  is 1.

generalize//@
      Solve[{Cos[x] == 2,  y/(x - 2*Pi - ArcCos[2]) == 3},  {x, y}]

Elimination of variables

A fairly common problem in computational math is to eliminate some variables from a set of equations. A frequently effective method is to use GroebnerBasis with a special monomial ordering. A common application is to the problem of finding implicit equations that define an algebraic variety given parametrically. For example, suppose we have the surface, given
parametrically by s and t, as

surf = {(((s^2) - 1 - (t^2))/((s^2) + (t^2) + 1)), 
        ((2 s)/((s^2) + (t^2) + 1)), ((2 s t)/((s^2) + (t^2) + 1))};

(You may recognize this as a common parametrization of a sphere). We clear denominators and add an auxiliary equation to prevent the 
denominator from vanishing. We then form a Gröbner basis with an ordering 
that eliminates the parameters and the auxiliary variable.

polys = {x ((s^2) + (t^2) + 1) - ((s^2) - 1 - (t^2)), 
        y ((s^2) + (t^2) + 1) - 2 s, z ((s^2) + (t^2) + 1) - 2 s t, 
        w ((s^2) + (t^2) + 1) - 1}GroebnerBasis[polys, {x, y, z}, {s, t, w}, 
      MonomialOrder -> EliminationOrder]

Here is a problem that is significantly more strenuous. It came from a Usenet news group question, from sci.math.symbolic, about a year ago. I was only able to handle it by finding an approximate Gröbner basis.

polys = {
        x - (110 (t^2) - 495 (t^3) + 1320 (t^4) - 2772 (t^5) + 5082 (t^6) - 
              7590 (t^7) + 8085 (t^8) - 5555 (t^9) + 2189 (t^10) - 374 (t^11))
            , y - (
            22 t - 110 (t^2) + 330 (t^3) - 1848 (t^5) + 3696 (t^6) - 
              3300 (t^7) + 1650 (t^8) - 550 (t^9) + 88 (t^10) + 22 (t^11))
            };

Timing[
      Short[poly250 = 
          First[GroebnerBasis[polys, {x, y}, t, 
              MonomialOrder -> EliminationOrder, 
              CoefficientDomain -> InexactNumbers[250]]]]]

Radical solutions vs. Root objects. vs. trigs for cubic equations

We are occasionally asked why Solve for cubic equations with real coefficients  returns results that contain I (that is, Sqrt[-1]) when all solutions for a given equation are real-valued. The short answer is that if the cubic is irreducible (does not factor) over the rationals then there is a theorem to the effect that two of the radical solutions will explicitly contain I; if all solutions are real these will manage to exactly cancel. All the same it is often nice to get solutions that do not have radicals at all. Below is an example.

rutes1 = Roots[(x^3) - 11 x + 2 == 0, x]

N[rutes1]

One can instead use algebraic numbers.

rutes2 = Roots[(x^3) - 11 x + 2 == 0, x, Cubics -> False]

N[rutes2]

Notice the ordering is not the same (no reason it ought to be, but sometimes we are asked about this).

You can even express the solution in terms of trigs and arctrigs. If you are into that sort of thing. It seems to work for the cubics I have tried; I am not sure if it works for cubics I have not tried (read: "caveat emptor"). Be warned that it is not a pretty sight.

rutes3  =  
      Apply[List,  Map[Last, rutes1 ]]  //  ComplexExpand

N[rutes3]

It is encouraging to note that Simplify can restore some flavor of radical result to this mess.

Simplify[Map[TrigToExp, rutes3]]

It is further encouraging to note that we can go from radicals to Root[...] results.

rutes1  /.  
      a_ == b_  :>  a  ==  RootReduce[b]

Finding approximate roots to polynomial equations

At the heart of Solve lies Roots. For numeric root finding there is the (undocumented) NRoots which in fact is called from Roots and elsewhere when needed. It uses the Jenkins-Traub algorithm. At times it can be a bit flaky. Whenever I suspect that it is misbehaving I will compare to the result of extracting eigenvalues of the polynomial's companion matrix. Below is code to do just this.

companionMatrix[poly_?PolynomialQ] := 
      Module[{vars, coeffs, newcoeffs, n, m}, 
        vars = Variables[poly]; 
If[
          ! ListQ[vars] || Length[vars] ≠ 1, 
          Return[]]; 
coeffs = 
          CoefficientList[poly, vars]; 
newcoeffs = 
          Drop[coeffs/(-Last[coeffs]), -1]; 
n 
          = Exponent[poly, First[vars]] - 1; 
m = 
          Prepend[IdentityMatrix[n], 
            Table[0, {n}]]; 
Transpose[
          Append[Transpose[m], newcoeffs]]]

Here is an example.

poly =  .7 + 3.2 x -  .3 (x^2) +  .9 (x^4) - 7.1 (x^5) + 1.3 (x^6) - 
        (x^7); MatrixForm[mat = companionMatrix[poly]]

Sort[Eigenvalues[companionMatrix[poly]]]

Here is the comparison with the built-in polynomial root-finder.

Sort[x /. Solve[poly == 0, x]]

Solving the cubic ...

Well, perhaps you've seen our excellent poster, "Solving the Quintic." The cubic is a classical result of the early sixteenth century (maybe fifteenth, my rustory is histy). Here we show how it is solved by radicals. Using Mathematica to simplify computational effort we will derive the general formula, the one you are told exists, but never wanted to look up in a reference. This method of derivation was explained to me by John D'Angelo. We begin with a generic cubic polynomial in x.

cubic = (x^3) + a (x^2) + b x + c;

We can get rid of the (x^2) term by a simple change of coordinates.When we 
finish, we then compensate by adding a/3 to each solution.

Collect[cubic /. x -> x - a/3 // Expand, x]

It thus suffices to find a formula for the cubic polynomial (x^3)+p*x+q. If 
you like (I do not), you can use what you saw above (that p is b - (a^2)/2 and q is c - a*b/3 + 2*(a^3)/27) to rewrite the final result in terms of a, b, and 
c rather than p and q.

newcubic = (x^3) + p x + q;

We now form a rational expression in t by the substitution x->(t+k/t).We will 
then clear denominators to get a degree six polynomial in t. After that we 
choose a "good" value for k, one that makes our polynomial essentially 
quadratic.

rat = Expand[newcubic /. x -> t + (k/t)]deg6poly = Collect[Expand[rat (t^3)], t]

We get the coefficient list,note that the degree four term has a coefficient 
linear in the new variable k, set it to zero and solve for k. We then check 
that this value also makes the degree two term vanish.We thus have a quadratic 
polynomial in (t^3).

coeffs = CoefficientList[deg6poly, t]newk = 
      k /. First[
          Solve[coeffs[[5]] == 0, k]]
          quad = 
      Expand[deg6poly /. k -> newk] /. 
        t -> (y^(1/3))

Now take a solution for y, using the quadratic formula. This gives rise to 
three solutions for t  (take any cube root of it, and then multiply that one 
by either  (((-1))^(2/3))or (((-1))^(4/3))).

ysol = 
      y /. Solve[quad == 0, 
            y][[1]]tsols = 
      {(ysol^(1/3)), (((-1))^(2/3)) (ysol^(1/3)), 
        (((-1))^(4/3)) (ysol^(1/3))}

We can now solve for x in terms of t.

xsols = (t + k/t) /. {k -> newk, t -> tsols}

newcubic /. x -> xsols // Together

Note that the solution will fail when one plugs in particular values for p 
and q that cause t to be zero (because we divide by t). This happens exactly 
when y is zero, which in turn happens when p is zero and Sqrt[q^2]==-q. This 
situation can only arise when polynomial has the (trivial) form (x^3)-q.

...and the quartic

If we begin with a monic quartic polynomial that has cubic term j (x^3) we can make the cubic term vanish by the change of variables x->(x-j/4). Thus we will simply begin with

quartic  =  x^4 + p*x^2 + q*x + r;

We next write it as a product of symbolic quadratics, that is, as ((x^2) + a x + b)*((x^2) - a x + c). If we can find {a,b,c} in terms of {p,q,r} then we can extract all four roots simply by two applications of the quadratic formula. We will first do the reverse, solving for {p,q,r} interms of {a,b,c}.

coeffs = 
        CoefficientList[
          x^4 + p*x^2 + q*x + r - (x^2 + a*x + b)*(x^2 - a*x + c), 
          x];
soln1 = Solve[coeffs == 0, {p, q, r}][[1]]

By inspection we observe that if we add (a^2) to p,and divide q by a, we get c+b and c-b.

Expand[{p + a^2, q/a} /. soln1]

Moreover,the difference of the squares of these is 4 b c, which is r. In other words, neweqn (defined below) is zero once the values for {p,q,r} are plugged into it.

neweqn = (p + a^2)^2 - (q/a)^2 - 4*r
Expand[neweqn /. soln1]

We clear the denominator, expand, and replace (a^2) by a new variable, asqr.
This gives a  cubic polynomial in asqr, with coefficients that are polynomials
in {p,q,r}. We solve it (we will use Solve, but one could use the cubic
solver outlined above). Note that one is free to choose any of the three
solutions for asqr (and also free to take either square root thereof down
below), so long as the choice is followed consistently.

cubic = 
      Expand[a^2*neweqn] /. 
        a^j_ -> asqr^(j/2)
asqrsol = 
      asqr /. Solve[cubic == 0, asqr][[1]]

We can now solve for {a,b,c} in terms of {p,q,r}.

{a, b, c} = {Sqrt[asqrsol], 1/2*(p + asqrsol - q/Sqrt[asqrsol]), 
          1/2*(p + asqrsol + q/Sqrt[asqrsol])};

We then extract our four solutions for x.

disc1 = Sqrt[a^2 - 4*b];
disc2 = Sqrt[a^2 - 4*c];
solns = 
      1/2*{-a + disc1, -a - disc1, a + disc2, a - disc2}

We will check these solutions by comparing with Solve. We also compare LeafCount values.

soln2 = 
        x /. Solve[x^4 + p*x^2 + q*x + r == 0, 
            x];
LeafCount[solns]
LeafCount[soln2]

We compare numerical results for these and two other methods.

solns /. {p -> 3, q -> -5, r -> 11} // N // 
      Sort
soln2 /. {p -> 3, q -> -5, r -> 11} // N // 
      Sort
x /. Solve[x^4 + 3*x^2 - 5*x + 11 == 0, x] // N // 
      Sort
x /. {ToRules[NRoots[x^4 + 3*x^2 - 5*x + 11 == 0, x]]} // 
      Sort

Note that the method will fail when the parameter a is zero. This happens
precisely when our quartic factors as ((x^2)  +  b)*((x^2)  +  c). In this case we
have a bi-quadratic equation, that is, an equation that is quadratic in (x^2).  Of course in this instance there is no need to use the heavier machinery above to extract the roots.

A tricky Solve example

This problem was suggested to me by Michael Trott. Given the infinite exponential "ladder" x^(x^(x^...)) == k, solve for x in terms of k. Noticing that the exponent of the first x is also just the infinite ladder we can rewrite this as the system {y == (x^y), y==k} (or, mentally substituting for y, we get (x^k)==k). We will solve the first system for x, eliminating y. It is a perversity of Solve that it can sometimes do this sort of thing; elimination of variables in a transcendental setting is an oddity. perhaps more unusual is that we must eliminate y; if we try to solve for both x and y we get a transcendental dependency.

Solve[{y == (x^y), y == k}, x, y]

Let us try this out for explicit values of k.

sol = 
      x /. First[Solve[{y == (x^y), y == E}, x]]

We will check this the old-fashioned way.

nsol = N[sol, 111]pow = nsol; Do[pow = (nsol^pow), {1000}]pow

pow  -  E

Not too bad. Some experimentation will indicate logarithmic convergence (iterate ten times as far, get one more digit correct). Now we try it for Pi.

sol = 
      x /. First[
          Solve[{y == (x^y), y == π}, x]]

sol = 
      x /. First[
          Solve[{y == (x^y), y == π}, x]]
nsol = N[sol, 111]
pow = nsol; Do[pow = (nsol^pow), {1000}]
pow - π

This time not so good. It turns out that further iteration will not help. So what went wrong? The problem lies in our assumption that the system has a solution. If it does not, then our "equivalent" system is really not equivalent at all.

More with LUDecomposition and linear equations

LUDecomposition of a square matrix is well used in the numeric linear algebra community. Perhaps less well known is that it can be applied in a modular setting. Below I show an implementation of a linear equation solver over the rationals that works as follows. We first renormalize to get a related system with integral solution. This is done simply by multiplying the right-hand-side by the determinant of the matrix. We will solve this new system modulo some prime, then iteratively "lift" the solution to be valid modulo increasing powers of that prime. In other words, we find a p-adic solution. As we go, we "correct" by symmetrizing with respect to these powers. Eventually we arrive at a solution that is correct over the integers. We then divide by the determinant to arrive at the solution for the original system.
The method is probabilistic in that we must choose a prime for which the matrix is non-singular. This will be true for all but at most finitely many primes (assuming it is nonsingular over the rationals). Note that
Ceiling[Log[N[mod], 2*matrixNorm[mat]*vectorNorm[rhs]]]
is an upper bound on the number of iterations through the main loop. In practice it will typically be less.

vectorNorm[vec_] := Apply[Plus, Map[Abs, vec]]
symmetrize[vec_, power_] := 
      Module[{len = Length[vec], halfpow = N[power/2]}, 
        Table[If[vec[[j]] > halfpow, vec[[j]] - power, 
            vec[[j]]], {j, len}]]
pAdicSolve[mat_?MatrixQ, rhs_?VectorQ] := 
      Module[{len = Length[mat], b, mod = Prime[2222], det, lud, sol, corr, 
          power}, lud = 
          LUDecomposition[mat, Modulus -> mod]; 
det = 
          Det[mat]; 
b = 
          det*rhs; 
power = 1; 
sol = 
          Table[0, {len}]; 
While[vectorNorm[b] > 0, 
          corr = symmetrize[LUBackSubstitution[lud, b, Modulus -> mod], 
              mod]; 
b = 
            1/mod*(b - mat . corr); 
sol += 
            power*corr; 
power *= 
            mod;
]; 
(1/det)*
          sol]

Here is an example.

mat = 
        Table[Random[
            Integer, {-(10^12), (10^12)}], {10}, {10}];
mod = Prime[2222];
b = 
        Table[Random[Integer, {-(10^12), (10^12)}], {10}];

Timing[ls1 = pAdicSolve[mat, b];]
Timing[ls2 = LinearSolve[mat, b];]
ls1 === ls2

mat = 
        Table[Random[
            Integer, {-(10^12), (10^12)}], {100}, {100}];
mod = Prime[2222];
b = 
        Table[Random[Integer, {-(10^12), (10^12)}], {100}];

Timing[ls1 = pAdicSolve[mat, b];]
Timing[ls2 = NewLinearSolve[mat, b];]
ls1 === ls2

Timing[ls1 = pAdicSolve[mat, b];]
Timing[ls2 = LinearSolve[mat, b];]
ls1 === ls2

LinearSolve[mat_, b_] :=  NewLinearSolve[mat, b]