4
$\begingroup$

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}]
$\endgroup$
8
  • 2
    $\begingroup$ I presume you know the Henderson-Hasselbalch equation? $\endgroup$ Commented Aug 24, 2012 at 13:57
  • 1
    $\begingroup$ 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. $\endgroup$
    – midtiby
    Commented Aug 24, 2012 at 13:58
  • $\begingroup$ You're using molars, I gather? $\endgroup$ Commented Aug 24, 2012 at 14:03
  • $\begingroup$ 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]]... $\endgroup$ Commented Aug 24, 2012 at 14:14
  • $\begingroup$ Can get something reasonable by plugging in a positive value for cHp, e.g. N@NSolve[(concentrations && equilibrium) /. {c2 -> x, cHp -> 2}, WorkingPrecision -> 50] $\endgroup$ Commented Aug 24, 2012 at 15:07

1 Answer 1

4
$\begingroup$

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}]

enter image description here

$\endgroup$
3
  • $\begingroup$ I don't think Solve makes use of $Assumptions which means that Assuming doesn't have any influence on Solve. $\endgroup$
    – Heike
    Commented Aug 24, 2012 at 20:15
  • $\begingroup$ @Heike, thanks. I'll edit my answer. $\endgroup$ Commented Aug 24, 2012 at 20:23
  • $\begingroup$ 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}]. $\endgroup$ Commented Aug 24, 2012 at 23:19

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.