Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 3 down vote favorite

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}]
share|improve this question
2  
I presume you know the Henderson-Hasselbalch equation? – J. M.♦ Aug 24 '12 at 13:57
1  
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? – J. M.♦ 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]]... – J. M.♦ 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
show 3 more comments

1 Answer

up vote 3 down vote accepted

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

share|improve this answer
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}]. – J. M.♦ Aug 24 '12 at 23:19

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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