# pH calculations, speeding up code involving Solve

I'm trying to calculate pH values of solutions containing water and acetic acid. The relevant equations have been defined and I now want to find a solution using Mathematica.

For a specific acid concentration (0.0001) I can do the following.

Solve[(concentrations && equilibrium && requirements) /. {c2 -> 0.0001}]


But I would prefer to solve the system of equations once and then insert the acid concentration. But the following code does not work

Solve[(concentrations && equilibrium && requirements)] /. {c2 -> 0.0001}


Which throws the following error

Solve::ratnz: Solve was unable to solve the system with inexact
coefficients. The answer was obtained by solving a corresponding
exact system and numericizing the result. >>"


My working code so far is shown here

concentrations = cHp == c1 + c3 && cOHm == c1 && cCH3COOm == c3 && cCH3COOH == c2 - c3;
equilibrium = 10^-14 == cHp cOHm && 1.1614 10^-5 == cCH3COOm cHp/cCH3COOH;
requirements = cHp > 0 && cOHm > 0 && cCH3COOH > 0 && cCH3COOm > 0;
Plot[-Log[10, cHp] /.
Solve[(concentrations && equilibrium && requirements) /. {c2 ->
x}], {x, 0, 0.00010}, PlotRange -> {2, 9}]

-
I presume you know the Henderson-Hasselbalch equation? –  Ｊ. Ｍ. Aug 24 '12 at 13:57
In that form I know it. But this example is meant as an example of what can be achieved with Mathematica. And as I can do the calculations by hand I expect that Mathematica can also do it if I ask it the right way. –  midtiby Aug 24 '12 at 13:58
You're using molars, I gather? –  Ｊ. Ｍ. Aug 24 '12 at 14:03
By the way, where'd you get your $K_a$ value for acetic acid? I get $\approx 1.74\times 10^{-5}$. Anyway, I can't see how you can do better than With[{c = 0.1}, -Log10[x] /. First@Solve[10^-4.76 == x^2/(c - x) && x > 0, x]]... –  Ｊ. Ｍ. Aug 24 '12 at 14:14
Can get something reasonable by plugging in a positive value for cHp, e.g. N@NSolve[(concentrations && equilibrium) /. {c2 -> x, cHp -> 2}, WorkingPrecision -> 50] –  Daniel Lichtblau Aug 24 '12 at 15:07

The following works for me, removing the requirements from Solve and filtering the solutions afterwards.

soln = Solve[concentrations && equilibrium, {cHp, cCH3COOm, cCH3COOH, c1, c3, cOHm}];


This gives three solutions, so select the one for which the requirements are satisfied at some reasonable value of c2.

soln = Select[soln, requirements /. Chop[# /. c2 -> 10^-5] &];


The solution can now be plotted:

Plot[-Log[10, cHp /. soln], {c2, 0, 0.00010}, PlotRange -> {2, 9}]


-
I don't think Solve makes use of \$Assumptions which means that Assuming doesn't have any influence on Solve. –  Heike Aug 24 '12 at 20:15
@Heike, thanks. I'll edit my answer. –  Simon Woods Aug 24 '12 at 20:23
Still, as I already mentioned, doing some of the work for Solve[] helps a lot: Plot[-Log10[x] /. First@Solve[10^-4.76 == x^2/(c - x) && x > 0, x], {c, 0, 10^-4}, PlotRange -> {2, 9}]. –  Ｊ. Ｍ. Aug 24 '12 at 23:19