I have a coupled chemical equilibria system consisting of 12 equations in 12 unknowns. I have managed use the Solve function to write a command that will find a particular solution at any fixed value of pH. I need to plot the value of each unknown across a given pH range. How can I do this?

  • $\begingroup$ Without seeing the code, it's hard to say (you can post the code here, and if it's large, try to condense it to just the most important part). $\endgroup$ – Jason B. Sep 8 '16 at 21:14

You sure have to wait sometimes to get a chemistry question around here lol. Perhaps this helps you in your PChem studies.

In the absence of your actual equations, let's take a similar example to the one you describe. Perhaps we know the initial concentration of phosphoric acid added , and we want to know the concentrations of all the various species in solution.

The chemical reactions involved are

enter image description here

That's four equations and six unknowns, but you can get two more equations by conservation of phosphate species, and conservation of protons (assuming an initial concentration of $10^{-7}$. You can feed all that to Solve, and if you also give conditions that the concentrations be positive then you get back only one solution. You've done that, but how do you plot it?

When you use Solve, it returns a list of replacement rules, like this

Solve[{x + y == 2, x - y == 2}, {x, y}]
(* {{x -> 2, y -> 0}} *)

To get the values, you need to use ReplaceAll, or it's infix form /.

{x, y} /. Solve[{x + y == 2, x - y == 2}, {x, y}]
(* {{2, 0}} *)

So here is the phosphoric acid concentration calculator:

  initConc_] :=
 {h, oh, h3PO4, h2PO4, hPO4, pO4} /. Quiet@Solve[
    {(h h2PO4)/h3PO4 == kA1,
     (h hPO4)/h2PO4 == kA2,
     (h pO4)/hPO4 == kA3,
     h oh == kW,
     h3PO4 + h2PO4 + hPO4 + pO4 == initConc,
     h == 3 pO4 + 2 hPO4 + h2PO4 + 10^-7,
     h > 0,
     h2PO4 > 0,
     h3PO4 > 0,
     pO4 > 0},
    {h, oh, h3PO4, h2PO4, hPO4, pO4}]

(as a PChem student, you'll probably want to include activity coefficients, but I'll leave them out for brevity). Let's test it by taking getting the pH for a 0.1M solution

(* {{0.0235359, 4.24882*10^-13, 0.0764642, 0.0235357, 
  6.20863*10^-8, 1.18923*10^-18}} *)
(* {{1.62827, 12.3717, 1.11654, 1.62827, 7.207, 17.9247}} *)

You can compare this to what you get from Alpha

WolframAlpha["pH 0.1M phosphoric acid", "Result"]
(* 1.63 *)

Looks good,* so how to plot it? I would collect the results in a list, and use ListLinePlot,

result = Transpose@Table[
   First@phosphoricAcidConcentrations[10^(-x)], {x, 0, 10, .1}];

 Frame -> True, FrameLabel -> {"-log(initial conc.)"},
 PlotLegends -> {"-log[H]", "-log[OH]", 
   "-log[\!\(\*SubscriptBox[\(H\), \(3\)]\)\!\(\*SubscriptBox[\(PO\), \
   "-log[\!\(\*SubscriptBox[\(H\), \
\(2\)]\)\!\(\*SubsuperscriptBox[\(PO\), \(4\), \(-\)]\)]", 
   "-log[\!\(\*SubsuperscriptBox[\(HPO\), \(4\), \(\(2\)\(-\)\)]\)]", 
   "-log[\!\(\*SubsuperscriptBox[\(PO\), \(4\), \(\(3\)\(-\)\)]\)]"}]

Mathematica graphics

Excuse the subscript-box gobbledegook.

*: But don't try it with very low concentrations, Alpha doesn't treat this as a polyprotic weak acid I think.

|improve this answer|||||

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.