Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

May I ask what is the best way to evaluate

NSolve[{a^10 E^-a - b^10 E^-b == 
   0, (362880 + 
       a (362880 + 
          a (181440 + 
             a (60480 + 
                a (15120 + 
                   a (3024 + 
                    a (504 + 
                    a (72 + a (9 + a))))))))) E^-a + (-362880 - 
       b (362880 + 
          b (181440 + 
             b (60480 + 
                b (15120 + 
                   b (3024 + 
                    b (504 + b (72 + b (9 + b))))))))) E^-b + 0.99 == 
   0}, {a, b}]

I tried several methods but my machine refuse to produce anything. Is the computational complexity really this much?

share|improve this question
By the way, though you've probably already done this, those numbers seem to match sequence A008279, a.k.a. permutation coefficients, a.k.a. falling factorials triangle. –  Andrew Cheong Apr 5 '14 at 6:09

2 Answers 2

up vote 7 down vote accepted

You can reformulate this as a objective function minimisation problem, where you minimise the sum of the squares of the left hand sides of your equations. Here is how to do this with NMinimize, trying all of the optimisation methods that are mentioned in the documentation:

{#, NMinimize[(a^10 E^-a - 
   b^10 E^-b)^2 + ((362880. + 
      a (362880. + 
         a (181440. + 
            a (60480. + 
               a (15120. + 
                a (3024. + 
                a (504. + 
                a (72. + a (9. + a))))))))) E^-a + (-362880. - 
      b (362880. + 
         b (181440. + 
            b (60480. + 
               b (15120. + 
                b (3024. + 
                b (504. + b (72. + b (9. + b))))))))) E^-b + 
   0.99)^2, {a, b}, Method -> #]} & /@
{"NelderMead", "DifferentialEvolution", "SimulatedAnnealing", "RandomSearch"}

which produces the output

{{"NelderMead", {0.965138, {a -> 1.96905, b -> 1.96881}}},
 {"DifferentialEvolution", {3.38813*10^-21, {a -> 1.69024, b -> -1.25858}}},
 {"SimulatedAnnealing", {1.11485*10^-20, {a -> 1.69024, b -> -1.25858}}},
 {"RandomSearch", {0.968532, {a -> 0.900447, b -> 0.774213}}}}

Both "DifferentialEvolution" and "SimulatedAnnealing" have found an accurate solution.

share|improve this answer
Ah, so that's how it's done! –  Andrew Cheong Apr 2 '14 at 14:35
I had another think about my solution when I read Andrew Cheong’s answer below. Define f1[a_,b_]:=<LHS #1>; f2[a_,b_]:=<LHS #2>, where "LHS" means "left hand side", then evaluate ContourPlot[{f1[a, b] == 0, f2[a, b] == 0}, {a, -50, 50}, {b, -50, 50}, ContourStyle -> {Red, Blue}, PlotPoints -> 100], where any solution corresponds to an intersection between the red and blue contours. Ignoring near-solutions where contours run very close to each other, there are 2 solutions: (1) close to (10,10) as noted by Andrew Cheong below, and (2) at {1.69024…, -1.25858…} as found in my answer above. –  Stephen Luttrell Apr 7 '14 at 18:58


There is another solution besides @StephenLuttrell's: $(10,9.99997819382)$.

At first I thought this was a mistake—namely a solution disguised as the trivial solution $(10,10)$ due to rounding issues—but it actually appears to be a real solution.


This answer has been updated in light of some glaring errors in my original post—namely in taking the sum of the logs of the squares instead of the log of the sum of the squares, and also in suggesting that zeroes would only appear as sharply ("infinitely") descending points, which was of course silly, since the lead-up to a discontinuity could take any shape, including "mild-looking" corners.

To rephrase, you're looking for the simultaneous zeroes of two functions, call them $f(a,b)$ and $g(a,b)$.

f[a_, b_] := a^10 E^-a - b^10 E^-b;
g[a_, b_] := (362880 + 
      a (362880 + 
         a (181440 + 
            a (60480 + 
               a (15120 + 
                  a (3024 + 
                    a (504 + 
                    a (72 + a (9 + a))))))))) E^-a + (-362880 - 
      b (362880 + 
         b (181440 + 
            b (60480 + 
               b (15120 + 
                  b (3024 + b (504 + b (72 + b (9 + b))))))))) E^-b + 

Here's a plot of $f$ and $g$ in blue and yellow, respectively:

enter image description here

For such unwieldy functions, it's easier to detect zeroes by squaring (or taking the absolute value of) the function, thereby removing all negative values, and apply minima-searching algorithms. And, to find the simultaneous zeroes of multiple functions, one simply minimizes the sum of their squares, plotted here:

enter image description here

Now, minima-searching algorithms don't care, but visually the above depiction isn't very useful to us, due to sheer scale. So, henceforth, we'll plot the Log of our functions (and sum of squares thereof) for the visual effect, but leave out the Log where applying algorithms. Here's the Log version of the above:

enter image description here

Alright—now we're talkin'! One of those spikes, marked by the green line, is the solution found by @StephenLuttrell. It's interesting to plot the sum of logs (as opposed to the log of sums) and zoom in on that region:

enter image description here

$f$ and $g$ create distinct "paths" of spikes, meeting at $(1.69024,-1.25858)$ as simultaneous zeroes.

By the way, plotting the log of squares (or absolute values) is useful in general for detecting zeroes visually, since $\log{x}$ as $x\to0$ tends toward $-\infty$, creating downward spikes at the zeroes.

Anyway, as @StephenLuttrell showed, NMinimize can be used to (attempt to) find the global minimum. However, can more solutions be found at other local minima? It seems we'd have to inspect each spike individually. Absent constraints, NMinimize searches for a global minimum, but making use of them, we can ask NMinimize to search, for example, only a square region around $(x_0,y_0)$, i.e. by specifying $$x_0-L<a<x_0+L\\y_0-L<b<y_0+L$$ where $L$ represents the size of each square (as the length of its sides).

To automate such inspections, we can "sweep" an area using two For loops, i.e. searching the area square by square:

enter image description here

As this process can take a while, as depicted above I created a DynamicModule that

  1. draws a red column indicating the area currently being searched,
  2. draws dotted red lines where (non-boundary) local minima were found, and
  3. draws solid green lines where $f^2+g^2<\text{threshold}$ for some $\text{threshold}$, e.g. $10^{-6}$.

See the bottom of this post for the source. Note that I started plotting f and g separately in order to take advantage of visual information—perhaps seeing intersecting paths of trenches can help us narrow down our search, visually.

Sweeping the area $(-2,-4)$ to $(6,4)$, searching in $0.2\times0.2$ squares, and applying all four NMinimize methods at each square, still only one point satisfied $f^2+g^2<10^{-6}$: the one already found by @StephenLuttrell. (The green line in the latest plot was not plotted manually—it was found and plotted by the search.)

Then, I multiplied everything by a factor of 100, i.e. now sweeping the area $(-200,-400)$ to $(600,400)$ and searching in $20\times20$ squares, though preserving the threshold at $10^{-6}$. This yielded a different solution, at approximately $(10,10)$!

enter image description here

Of course, it overrode @StephenLuttrell's solution, because it happened to be the global minimum in the $20\times20$ square that includes the $0.2\times0.2$ square containing @StephenLuttrell's solution. Manually, I found the more precise solution to be $(10, 9.99997819382)$:

f[10, 10]
(* 0 *)

g[10, 9.99997819382]
(* -1.19762*10^-7 *)

Epilogue (and not the Graphics kind...)

Just for fun, I zoomed in on $(10,10)$ as I had done with @StephenLuttrell's solution:

enter image description here

(The green line's actually a bit off, due to an apparent precision limit on Plot3D.)

What is this crazy thing! Might there be another solution nearby, where the other intersection is? Adjusting my code's parameters to search $(9.99995,9.99995)$ to $(10.00005,10.00005)$, searching in $.00001\times.00001$ squares, and now adjusting the threshold to $10^{-11}$, it was immediately apparent that there would be no other solutions; $g$ had only one trail of spikes—the "parallel" trail belonged to $f$:

enter image description here


   searchRadius = 20,
   searchOverlap = 1,
   searchXmin = -200,
   searchXmax = 600,
   searchYmin = -400,
   searchYmax = 400,
   plotZmin = -1000,
   plotZmax = 1000,
   currentSearchArea = {},
   localMinima = {},
   localMinimaBelowThreshold = {},
   threshold = 10^-6,
   algorithms = {"NelderMead", "DifferentialEvolution", 
     "SimulatedAnnealing", "RandomSearch"},
   i, j, imin, imax, jmin, jmax, x, y
     {Log[f[x, y]^2], Log[g[x, y]^2]},
     {x, searchXmin, searchXmax},
     {y, searchYmin, searchYmax},
     PlotRange -> {
       {searchXmin - searchRadius - searchOverlap/2, 
        searchXmax + searchRadius + searchOverlap/2},
       {searchYmin - searchRadius - searchOverlap/2, 
        searchYmax + searchRadius + searchOverlap/2},
       {plotZmin, plotZmax}},
     PlotStyle -> {
       {Opacity[0.5], Blue},
       {Opacity[0.5], Yellow}
     MeshStyle -> Gray,
     ImageSize -> Full
    Graphics3D[{Red, Opacity[0.7], EdgeForm[None], 
    Graphics3D[{Red, Dotted, 
        Line[{{t[[1]], t[[2]], plotZmin}, {t[[1]], t[[2]], 
           plotZmax}}], {t, localMinima}]]}],
        Line[{{t[[1]], t[[2]], plotZmin}, {t[[1]], t[[2]], 
           plotZmax}}], {t, localMinimaBelowThreshold}]]}]
  For[i = searchXmin, i <= searchXmax, i += 2*searchRadius,
   For[j = searchYmin, j <= searchYmax, j += 2*searchRadius,
    imin = i - searchRadius - searchOverlap/2;
    imax = i + searchRadius + searchOverlap/2;
    jmin = j - searchRadius - searchOverlap/2;
    jmax = j + searchRadius + searchOverlap/2;
    currentSearchArea = {{imin, jmin, plotZmin}, {imax, jmax, 
       {x, y} = {a, b} /. #;
       If[x != imin && x != imax && y != jmin && y != jmax,
        If[f[x, y]^2 + g[x, y]^2 < threshold,
          Print[{x, y, f[x, y], g[x, y]}];

          localMinimaBelowThreshold = 
           Append[localMinimaBelowThreshold, {x, y}],
          localMinima = Append[localMinima, {x, y}]
        ]) &,
        NMinimize[{f[a, b]^2 + g[a, b]^2, {imin < a < imax, 
            jmin < b < jmax}}, {a, b}, Method -> #][[2]] & /@ 
  currentSearchArea = {{0, 0, 0}, {0, 0, 0}};
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.