1
$\begingroup$

Dear Mathematica experts. I am a chemist looking to use Mathematica to solve a system of equilibrium equations. I have a system of 8 polynomial (I think) equations with 8 unknowns, yet when I use nsolve I end up with an empty set {}. Does anyone have any advice? Thank you in advance for your help. ;)

NSolve[{x1*x7*10^6 == x3, x1*x8*1.67*10^8 == x5, x3*x8*10^6 == x4, 
  x5*x8*10^6 == x6, x1*x8*10^6 == x2, 
  10 == x1 + x3 + x4 + x5 + x2 + x6, 500 == x3 + x4 + x7, 
  12 == x5 + x2 + x4 + 2*x6 + x8}, {x1, x2, x3, x4, x5, x6, x7, 
  x8}, Reals]
$\endgroup$
  • $\begingroup$ The code work for me. What is your version of Mathematica? $\endgroup$ – Kattern May 28 '15 at 11:43
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey May 28 '15 at 11:51
  • $\begingroup$ kattern, I am using the trial version of Mathematica. Thank you for your help. $\endgroup$ – Katy May 28 '15 at 12:20
  • 1
    $\begingroup$ I suspect NSolve at WorkingPrecision->Automatic is getting smallish residuals from results that are not so great (possibly a consequence of scaling differences in the input). If you do say NSolve[{...,WorkingPrecision -> 100] you may get a result more to your liking. For this to be quietly effective you will need exact (or at least high precision) input, which can be done simply by replacing the coefficient 1.67 with 167/100 (or can do SetPrecision[...,Infinity] on the input system). $\endgroup$ – Daniel Lichtblau May 28 '15 at 15:16
1
$\begingroup$

If those $x_i$ values represent concentrations of chemical species, then of course they don't just have to be real, they also have to be positive numbers or very very close to zero, accounting for numerical precision errors.

You can impose that restriction in Solve directly. Solve will complain trying to solve your problem with inexact coefficient (i.e. 1.68*10^8), so you can artificially make that number exact as 168*10^6 in the following equations.

In short:

solutions = Solve[
   {
    x1*x7*10^6 == x3, x1*x8*167*10^6 == x5, x3*x8*10^6 == x4,
    x5*x8*10^6 == x6, x1*x8*10^6 == x2, 10 == x1 + x3 + x4 + x5 + x2 + x6, 
    500 == x3 + x4 + x7, 12 == x5 + x2 + x4 + 2*x6 + x8,
    x1 >= 0, x2 >= 0, x3 >= 0, x4 >= 0, x5 >= 0, 
    x6 >= 0, x7 >= 0, x8 >= 0
   }, 
   {x1, x2, x3, x4, x5, x6, x7, x8}
];

These expressions are quite ugly, but you probably don't care about their mathematical form, so we can just proceed to obtain a numerical value using N:

N[{x1, x2, x3, x4, x5, x6, x7, x8} /. solutions]

{{3.37165*10^-14, 1.73224*10^-8, 0.0000165712, 8.51373,
  2.89284*10^-6, 1.48625, 491.486, 0.513767}}
$\endgroup$
  • $\begingroup$ Hi Marco, this looks great as well. i'm going to try it out now. Thank you so much! :) K $\endgroup$ – Katy May 28 '15 at 13:15
  • $\begingroup$ Sorry I have a silly question. Is it possible to generate those 8 solutions in some kind of row or column table so I can just paste it into excel? Right now I have to use my mouse to highlight and copy and paste each number into a cell separately. Thanks a bunch! $\endgroup$ – Katy May 28 '15 at 13:23
  • 1
    $\begingroup$ @Katy try this: Add Column@(CForm /@ First@%) after the expression that evaluates what you want to copy (e.g. the N expression in my answer above). This will give a column-formatted list of values in scientific "E" notation that Excel will understand. Select all the values at once, right click and select "Copy As..." -> "Plain Text", then paste in Excel. It should give you a nice column of values in the spreadsheet. For longer lists, look into generating a CSV or XLSX file directly using the Export command: in that case the formatting would be adjusted automatically. $\endgroup$ – MarcoB May 28 '15 at 13:50
  • $\begingroup$ MarcoB, you are a genius! Thank you!!! $\endgroup$ – Katy May 28 '15 at 14:07
1
$\begingroup$
$Version

"10.1.0 for Mac OS X x86 (64-bit) (March 24, 2015)"

eqns = {x1*x7*10^6 == x3, x1*x8*1.67*10^8 == x5, x3*x8*10^6 == x4, 
    x5*x8*10^6 == x6, x1*x8*10^6 == x2, 10 == x1 + x3 + x4 + x5 + x2 + x6, 
    500 == x3 + x4 + x7, 12 == x5 + x2 + x4 + 2*x6 + x8} // Rationalize;

var = Variables[Level[eqns, {-1}]];

As you indicated, NSolve does not return a solution for the Reals

NSolve[eqns, var, Reals]

{}

However, you can use FindInstance

soln1 = FindInstance[And @@ eqns, var, Reals, 5] // N;

Verify solutions

test = And @@ eqns /. soln1;

Select valid solutions

soln2 = Pick[soln1, test]

{{x1 -> -502.073, x2 -> 0.0718752, x3 -> 500.072, x4 -> -0.0715886,
x5 -> 12.0032, x6 -> -0.00171833, x7 -> -9.96013*10^-7, x8 -> -1.43157*10^-10}, {x1 -> -4.0817*10^-9, x2 -> 2.44899*10^-8, x3 -> -2.00003, x4 -> 12.0001, x5 -> 4.08981*10^-6, x6 -> -0.0000245386, x7 -> 490., x8 -> -5.99993*10^-6}, {x1 -> 6.08952*10^-16, x2 -> -7.026*10^-9, x3 -> 3.0663*10^-7, x4 -> -3.53786, x5 -> -1.17334*10^-6, x6 -> 13.5379, x7 -> 503.538, x8 -> -11.5379}, {x1 -> 3.37165*10^-14, x2 -> 1.73224*10^-8, x3 -> 0.0000165712, x4 -> 8.51373, x5 -> 2.89284*10^-6, x6 -> 1.48625, x7 -> 491.486, x8 -> 0.513767}}

$\endgroup$
  • $\begingroup$ Thank you Bob! Another option! that's why you guys are the maths experts and I'm just a chemist! :) $\endgroup$ – Katy May 28 '15 at 14:08
  • $\begingroup$ Don't jump to that conclusion too soon, @Katy, the author of the answer you accepted is also a chemist… $\endgroup$ – J. M. is in limbo May 28 '15 at 16:20
0
$\begingroup$

Using FindRoot with a wild guess of starting values does return a solution:

FindRoot[{
  x1*x7*10^6 == x3,
  x1*x8*1.67*10^8 == x5, 
  x3*x8*10^6 == x4, 
  x5*x8*10^6 == x6, 
  x1*x8*10^6 == x2,
  10 == x1 + x3 + x4 + x5 + x2 + x6,
  500 == x3 + x4 + x7, 
  12 == x5 + x2 + x4 + 2*x6 + x8
  }, {{x1, 0}, {x2, 0}, {x3, 0}, {x4, 0}, {x5, 0}, {x6, 0}, {x7, 0}, {x8, 0}}
]

returns

{x1 -> -4.0242791928197337`*^-16, x2 -> -5.0298429418384916`*^-9, 
 x3 -> 1.6400641951497724`*^-6, x4 -> 20.498740074631854`, 
 x5 -> -8.399837712870253`*^-7, x6 -> -10.498740869682432`, 
 x7 -> 479.50125828530395`, x8 -> 12.498742509746629`}

which may or may not be an acceptable root. (In particular, this makes x1*x7*10^6 == -1.92965*10^-7 but x3 == 1.64006*10^-6 for the first equation, which may or may not be an acceptable discrepancy.)

It's important to note that a system of $n$ independent polynomial equations in $n$ variables need not have any real solutions. (As a starting example, take $x^2+1=0$ in one variable.) If you have faith that your equations correctly describe the equilibrium of a system and you have good reasons to think the equilibrium exists, then you can be confident that the equations will have solutions. If Mathematica is having trouble finding them or you do not like its results, the first thing to try is to use your intuition about your chemical system to come up with a reasonable first guess for the solution, and then feed this to FindRoot.

For your particular system, one thing which seems to be confusing Mathematica is that your variables have wildly different magnitudes, which is evident by the presence of very large numbers in the equation. If you have definite expectations for the order of magnitude of each variable at the equilibrium point, try rescaling them so that all the variables will be of order 1.

Beyond that, if algebraic methods don't really work, you can try scaling back up to your full rate equations and then use NDSolve to make them relax to the equilibrium point.

$\endgroup$
  • $\begingroup$ Thank you episanty!! It worked! Hurray! Thanks for pointing out the large and small values. I actually made an embarrassing error with my original equation, the 10, 500 and 12 should all be multipled by 10^-6, but Nsolve didn't work. I took your suggestion and converted all units to micromolar scale, to give numbers ranging from 10 to 167, but still Nsolve returned nothing. However find root was successful! I wonder what is the functional difference between find root and Nsolve that allows find root to work but Nsolve not? thank you so much!!! you are a lifesaver! $\endgroup$ – Katy May 28 '15 at 12:25
  • $\begingroup$ @Katy, Is it acceptable for your formulation of the problem to have negative concentration values, such as the x6 found above by FindRoot? If not, you can either give FindRoot a better guess of starting points, if you have one, or still solve the system exactly (take a look at my answer below). $\endgroup$ – MarcoB May 28 '15 at 12:51
  • $\begingroup$ hi Marco, you are correct, we cannot have negative concentrations. I overcame this by giving it better guess of starting points. Thank you very much for your help. I will take a look at your answer below as well. $\endgroup$ – Katy May 28 '15 at 13:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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