Short version: How do I have Solve handle a system of equations where some of the equations may be conditional?

Simplified Problem:

I am solving the following system of equations:

  1. a * b = 10^-14
  2. b^2 * c = 5.02 * 10^-6
  3. b * d = 25


  1. a -> a0 + x1
  2. b -> b0 + x1 + 2 * x2 + x3
  3. c -> c0 + x2
  4. d -> d0 + x3

a0, b0, c0, and d0 are known constants. I am solving for x1, x2, and x3.

But... x2 has a maximum allowed value, and if that maximum value is reached, then equations 2 is thrown out. x3 also has a maximum value, and if that maximum value is reached, equation 3 is thrown out.

For testing, suppose that a0, b0, c0, and d0 are all zero. For maximum values of x2 and x3, let’s try the following scenarios:

A) x2max = 1, x3max = 50 (neither limit will kick in)

B) x2max = 1, x3max = 1 (x3 limit will kick in)

C) x2max = 10^-10, x3max = 50 (x2 limit will kick in)

D) x2max = 10^-10, x3max = 1 (both limits kick in)

How do I set up that kind of scenario in a Solve statement?

Context: Suppose that I have three chemical equilibria:

enter image description here

These are governed by equilibrium constants:

enter image description here

Here, in chemistry-standard notation, the square brackets mean "concentration of", and the M is a unit of concentration (moles per liter).

A typical strategy would have me identify candidate solutions to the problem by taking initial, known concentrations of all of the species marked (aq), and adjusting them based on the "extent of reaction" for each of the equilibria. So each molecule of water that undergoes the top reaction would decrease the amount of water by one molecule, and increase the amount of the two products each by one molecule. So if we label the extent of reaction for the three reactions x1, x2, and x3, the final concentration of OH-(aq) would be the initial concentration plus x1 + 2 x2 + x3. I plug similar expressions in to the three equilibrium constant expressions and solve for x1, x2, and x3. I would do something like:

Solve[{H OH == 10^-14,OH^2 Ca==5.02*10^-6,OH Na==25}/.
        {H->Hinit+x1, OH->OHinit+x1+2 x2+x3, 

I can then plug the results back into the concentration expressions and pick the only option for which all of the concentrations are positive.

But suppose that I only have a limited amount of one of the solids, say NaOH. That means that the last equation doesn't hold, because all of the solid dissolves before reaching the massive concentrations that would be necessary for that equilibrium to be in-play. This is kind of like a boundary condition. If I am doing this as a one-off problem, it is simple enough to calculate the concentration that you get when all that I do have dissolves, use that for my starting concentrations for the other equilibria, and then throw out the third equation since it doesn't hold. But I am instead trying to develop some general routines that can recognize this kind of situation and address it.

Back in my procedural programming days, I would probably have done exactly that. I would come up with a candidate solution, then tested to see if I ran out of one of the solids, then modify the starting conditions and the list of equations to solve accordingly, but I have seen so many magical ways for Mathematica to simplify such approaches that I feel there must be an elegant way to do it here. Suggestions?

Code Example:

Let's say I set up my equilibria as follows:

equilibria = {Equ[H2O, H + OH, Quantity[10^-14, "Molar"^2]],
              Equ[CaOH2, 2 OH + Ca, Quantity[5.02*10^-6, "Molar"^3]], 
              Equ[NaOH, OH + Na, Quantity[25, "Molar"^2]]};

Here, Equ is a function that I never define; it is just used to hold the structure. The first parameter is the reactant side; the second parameter is the product side; and the third parameter is the equilibrium constant. Next I define some functions I will be using to make my life easier, essentially extracting various aspects of equilibria:

Reactants[Equ[a_, b_, ___]] := Variables[a]
Products[Equ[a_, b_, ___]] := Variables[b]
ReactantStoichiometry[Equ[a_, b_, ___]] := Coefficient[a, Variables[a]]
ProductStoichiometry[Equ[a_, b_, ___]] := Coefficient[b, Variables[b]]
Coeff[Equ[a_, b_, ___], species_] := Coefficient[b, species] - 
                                     Coefficient[a, species]
Species[Equ[a_, b_, ___]] := Union[Variables[a], Variables[b]]
EqConst[Equ[a_, b_, k_]] := k

This gives me the reactant species, the product species, their coefficients (known as the "Stoichiometry") the coefficient of a specific species (positive if a product, negative if a reactant), a full list of the species, and the equlibrium constant.

I also need the states of matter (s = solid, l = liquid, aq = aqueous solute), because l gives me the volume (the denominator in the concentration calculation), s is the solids we might run out of, and aq are the species that are factors in the equilibrium constants:

State = <|H2O -> "l", H -> "aq", OH -> "aq", CaOH2 -> "s", 
          Ca -> "aq", NaOH -> "s", Na -> "aq"|>;

Now I can automatically determine the equilibrium expressions:

KeqConc[rxn_] :=
      Times @@ (((If[MemberQ[{"sol", "aq"}, State[#]], 
      conc[#], 1]) & /@ Products[rxn])^ProductStoichiometry[rxn])/
      Times @@ (((If[MemberQ[{"sol", "aq"}, State[#]], 
      conc[#], 1]) & /@ Reactants[rxn])^ReactantStoichiometry[rxn])

Here, the concentration of a species is represented by conc[species].

KeqConc[#] & /@ equilibria
(* {conc[H] conc[OH], conc[Ca] conc[OH]^2, conc[Na] conc[OH]} *)

And we see the equilibrium expressions are determined correctly.

We need a generic way to describe a sample. We'll do it as a list of components and their amounts. We'll take 0.1 liter of water and add 0.1 mole CaOH2 and 0.1 mole NaOH to it.

sample = {{H2O, Quantity[0.1, "Liters"]}, 
          {CaOH2, Quantity[0.1, "Moles"]}, 
          {NaOH, Quantity[0.1, "Moles"]}};

We can automatically determine what all of the solutes are--the ones designated "aq" or "sol" ("sol" is used when the solvent isn't water):

soluteList = Union[Cases[Union @@ (Species[#] & /@ equilibria), 
     _?(StringMatchQ[State[#], "sol" | "aq"] &)], 
  Cases[sample[[All, 1]], _?(StringMatchQ[State[#], "sol" | "aq"] &)]]
(* {Ca, H, Na, OH} *)

That may look unnecessarily complicated, but it is possible to have species in the samples that do not participate in any of the equilibria, so my code has to account for that.

I also have to calculate the volume:

volume = Total[Transpose[Select[sample, StringMatchQ[State[#[[1]]], "l"] &]][[2]]]
(*  0.1 L  *)

We also need to determine the initial conditions (before any of the equilibria have occurred). In this case, we set it up so that they are all zero, but that need not have been the case. Here is the general way to do it:

initConditions = (initialconc[#] -> ((Total[
         sample, {#, _}], {{#, 
          Quantity[0, "Moles"]}}]\[Transpose][[2]]])/volume)) & /@ soluteList
(*  {initialconc[Ca] -> Quantity[0., ("Moles")/("Liters")],
     initialconc[H] -> Quantity[0., ("Moles")/("Liters")], 
     initialconc[Na] -> Quantity[0., ("Moles")/("Liters")], 
     initialconc[OH] -> Quantity[0., ("Moles")/("Liters")]}  *)

We now determine the final concentrations based on the extent of reaction for each equilibrium:

finalConc = (#1 -> #2) & @@@ ({(conc[#] & /@ 
   soluteList), (initialconc[#] & /@ soluteList) + 
  Total[(Thread[Coeff[#, soluteList]] & /@ equilibria) (x /@ 

We have Solve find our candidate solutions:

candidateSolutions = 
  Quiet[Solve[(KeqConc[#] == EqConst[#]) & /@ equilibria /. 
         finalConc /. initConditions, 
       (x /@ Range[Length[equilibria]]), Reals]]

enter image description here

Only one of which has non-negative concentrations:

Cases[finalConc /. candidateSolutions /. initConditions, 
      _?(AllTrue[QuantityMagnitude[(conc[#] & /@ soluteList) /. #], 
     (# >= 0) &] &)]
(*  {{conc[Ca] -> Quantity[0.0002008, ("Moles")/("Meters")^3], 
      conc[H] -> Quantity[2.*10^-12, ("Moles")/("Meters")^3], 
      conc[Na] -> Quantity[5000., ("Moles")/("Meters")^3], 
      conc[OH] -> Quantity[5000., ("Moles")/("Meters")^3]}}  *)

But, as you see, this has not taken the limited amounts of the solids into account, which is the thrust of my question. I suppose one could set up the extent of reaction to have a maximum value somehow? I'm really not sure where to start.

For anyone that is curious, this is part of a package I am developing that looks like this:

enter image description here

  • 2
    $\begingroup$ I'm a bit confused by all this. It presumably represents a reaction network. What type kinetics (e.g. mass action)? I'd expect some ODEs somewhere, and finding equiliibria amounts to solving the algebraic system that results when the derivatives are set to zero. $\endgroup$ Commented Jun 26, 2019 at 15:33
  • $\begingroup$ No, it is not modeling kinetics. I have a separate set of routines for doing that. Equilibrium problems are a question of thermodynamics, not kinetics. No ODEs needed. Equilibrium says that the product of the concentrations of the products, divided by the product of the concentrations of the reactants is equal to a constant. $\endgroup$ Commented Jun 26, 2019 at 18:21
  • $\begingroup$ I have edited the question to provide a simplified version of the problem. $\endgroup$ Commented Jun 26, 2019 at 18:38
  • $\begingroup$ Can you provide example sets of values for a0, b0, c0, d0 and the maximum values of x2 and x3? $\endgroup$
    – Carl Woll
    Commented Jun 26, 2019 at 18:58
  • $\begingroup$ @CarlWoll: Info added to OP. $\endgroup$ Commented Jun 26, 2019 at 19:20

1 Answer 1


Here is a possibility. Your equations:

    a b == 10^-14,
    b^2 c == 5.02*^-6,
    b d == 25
} /. {a -> a0+x1, b -> b0 + x1 + 2 x2 + x3, c -> c0 + x2, d -> d0 + x3}

{(a0 + x1) (b0 + x1 + 2 x2 + x3) == 1/ 100000000000000, (c0 + x2) (b0 + x1 + 2 x2 + x3)^2 == 5.02*10^-6, (d0 + x3) (b0 + x1 + 2 x2 + x3) == 25}

Using Solve for your second case:

Block[{a0 = 0, b0 = 0, c0 = 0, d0 = 0, r},
    r = Solve[
            x2==1 || (x2<1 && eqns[[2]]),
            x3==1 || (x3<1 && eqns[[3]])
        {x1, x2, x3},
    MaximalBy[r, Count[eqns /. #, True]&]

{{x1 -> 9.9999*10^-15, x2 -> 5.0199*10^-6, x3 -> 1.}}

  • $\begingroup$ Something isn't working there, because all of those solutions fail the criteria of x2 and x3 being capped at 1. $\endgroup$ Commented Jun 26, 2019 at 20:31
  • $\begingroup$ We get a lot closer with this, which gives only two solutions, the second of which is correct: examplesolutions = Block[{a0 = 0, b0 = 0, c0 = 0, d0 = 0}, FindInstance[ And[eqns[[1]], x2 == 1 || (eqns[[2]] && x2 < 1), x3 == 1 || (eqns[[3]] && x3 < 1), x1 > 0, x2 > 0, x3 > 0], {x1, x2, x3}, Reals, 10]] // Column $\endgroup$ Commented Jun 26, 2019 at 20:36
  • $\begingroup$ @KevinAusman See update. $\endgroup$
    – Carl Woll
    Commented Jun 26, 2019 at 20:53
  • $\begingroup$ Oooo... I like your MaximalBy approach! Very nice! And much better than the approach I was pursuing. Thank you, this looks to be exactly what I need! $\endgroup$ Commented Jun 26, 2019 at 21:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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