The system is an aqueous solution of sulfuric acid and sulfur dioxide. I know how much sulfur there is and I know the pH. The system is in equilibrium. I would like to use Octave to solve the system.

I set up the following equations, ignoring $\ce{SO3^2-}$ because the pH is around 1.5.

By definition: $$ \mathrm{pH} = -\log [\ce{H+}]$$

Conservation of mass: $$[\ce{H2SO4}] + [\ce{H2SO3}] + [\ce{HSO4-}] + [\ce{SO4^2-}] + [\ce{HSO3-}] = M $$

Henderson-Hasselbalch equations for the system: \begin{align} \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{HSO3-}) + \log \frac{[\ce{HSO3-}]}{[\ce{H2SO3}]}\\ \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{SO4^2-}) + \log \frac{[\ce{SO4^2-}]}{[\ce{HSO4-}]}\\ \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{HSO4-}) + \log \frac{[\ce{HSO4-}]}{[\ce{H2SO4}]}\\ \end{align}

I then defined the solution vector to have the following order: $$ x = [ [\ce{H+}], [\ce{HSO4-}], [\ce{SO4^2-}], [\ce{HSO3-}], [\ce{H2SO4}], [\ce{H2SO3}] ]$$

and my code for Octave is attached. For some reason I always get negative values for some of the results.

function y = phCalc(x)
   pH = 1.87;
   pKa1 = -5;   % this is a guess to drive H2SO4 to ionization in this system
   pKa2 = 1.987;% from a textbook
   pKa3 = 1.857;% from a text book
   M = 0.0012 / 32 * 1000 ; % total sulfur moles per volume
   y(1) = 10^(-pH) - x(1);
   y(2) = x(5) * 10^(pH-pKa1) - x(2);
   y(3) = x(2) * 10^(pH-pKa2) - x(3);
   y(4) = x(1) - x(2) - x(3) - x(4);
   y(5) = x(2) + x(3) + x(4) + x(5) + x(6) - M;
   y(6) = x(4) / 10^(pH-pKa3) - x(6);
  • 2
    $\begingroup$ This question seems to be about programming and not about chemistry. Therefore, it might be considered off-topic here. $\endgroup$
    – chipbuster
    Commented Aug 10, 2015 at 22:26
  • 1
    $\begingroup$ If you can separate the chemistry question from the Octave debugging question, that will help us help you work towards solving your problem. $\endgroup$ Commented Aug 10, 2015 at 23:04
  • 3
    $\begingroup$ I don't think this question is about programming. Presenting the code is actually a nice unambiguous way to describe the chemical calculation that is being attempted. The problem seems to be with the chemistry rather than with the code... $\endgroup$
    – Curt F.
    Commented Aug 11, 2015 at 1:06
  • 1
    $\begingroup$ Welcome to Chemistry.SE! Take the tour to get familiar with this site. Mathematical formulae can be formatted using $\LaTeX$ syntax.|| I think this question is well on topic here, it is just like a homework question, where the attempt is writing out code. However, I might be able to understand the chemistry behind, but I cannot understand the code. I think it is quite necessary to explain all the variables and constants in use. A chemical equation of the equilibrium would also help. $\endgroup$ Commented Aug 11, 2015 at 9:22
  • 2
    $\begingroup$ I think I see the problem now. For this system to approach this pH, M has to float a bit more, and cannot be solved using that number. The lowest value for this M is around 2.5. At this pH, this system only has positive solutions between [0.00031..0.00085]*32/1000. Apparently there should be expressed another constraint, or maybe the problem is overconstrained. What do you think? $\endgroup$
    – JohnS
    Commented Aug 11, 2015 at 17:05

2 Answers 2


Your code seems to be calculating the differences between a guessed concentration vector x and a calculated concentration vector. These differences in concentrations between guessed and calculated are y. For example, the guessed pH is apparently x(1), and the calculated pH is 10^(-pH). The difference between them is y(1).

Thus, its no surprise that the y values are sometimes negative. To solve the problem you need to feed your function to a numerical solver such as fsolve. That function will take your guessed solution x, other parameters such as pH, the various pKas, and M, and your function phCalc as input, and return the value for x which makes all the y variables zero. That solution is the one you want.

On the chemistry side, you mention sulfur dioxide, but I don't see any equations for the solubility of sulfur dioxide gas in water. Additionally, sulfur dioxide is reduced relative to sulfuric acid, so there would need to be an additional reaction, and parameter (extent of reaction or degree of reduction or something) that describes the relative amounts of reduced sulfur (as sulfur dioxide, bisulfite, or sulfite) vs. oxidized sulfur (as sulfate, bisulfate, sulfuric acid, or sulfur trioxide).

There is also the issue of charge balance: what concentration of inert cations (e.g. sodium or potassium or whatever) are you assuming? Charge balance must hold in bulk aqueous solution. If you don't have any sodium / potassium etc., then you should include the self ionization of water as an equation, hydroxide as a species, and y(7) = x(1) - x(2) - 2*x(3) - x(4) - x(OH) as an equation. If you do have sodium/potassium/whatever at significant concentration, you can neglect hydroxide and just go with the equation y(7) = x(1) - x(2) - 2*x(3) - x(4) + x(Na) where x(Na) is the cation concentration, which I have assumed is sodium.


  • $\begingroup$ That's right. fsolve() gives an x to the user's function in an attempt to minimize the result (the resultant vector should be zero when a solution is found). That's just the way it's written. This code works correctly if I try to calculate only [H+], for example, by setting y(2) = x(2); y(3) =x(3); y(4)=x(4); y(5)=x(5); y(6)=x(6), and leaving y(1) as in the code above. So did I make a mistake in the chemistry? $\endgroup$
    – JohnS
    Commented Aug 11, 2015 at 13:55
  • $\begingroup$ There are a few things about the chemistry in your equations that I'm not sure of. Let me edit my answer to expand. $\endgroup$
    – Curt F.
    Commented Aug 11, 2015 at 15:48
  • $\begingroup$ Actually I don't have any sodium or potassium. The charge balance I forgot to mention but was included in the calculation of y(4). $\endgroup$
    – JohnS
    Commented Aug 11, 2015 at 17:46
  • $\begingroup$ Thanks for the help. I think what I did wrong was assume a mass / pH for this system which is nonrealistic. $\endgroup$
    – JohnS
    Commented Aug 11, 2015 at 17:47

What you want to solve, is the following system: $$ \frac{k_W}{\ce{[H+]}} - \ce{[H+]} + \frac{\ce{[H2SO4]0}\left(\ce{[H+]}~k_{a1}+2~k_{a1}k_{a2}\right)} {\ce{[H+]}^2+\ce{[H+]}~k_{a1}+k_{a1}k_{a2}} + \frac{\ce{[H2SO3]0}\left(\ce{[H+]}~k_{a3}+2~k_{a3}k_{a4}\right)} {\ce{[H+]}^2+\ce{[H+]}~k_{a3}+k_{a3}k_{a4}} = 0 $$ for both concentrations of $\ce{[H2SO4]0}$ and $\ce{[H2SO3]0}$ with the constraint that both concentrations add up to $\mathrm{0.0375~mol/L}$.

I have no idea how to solved this with Octave (see below), but as I am more used to Mathematica, that was my first approach to see where the error came from. The problem that you have with those negative solutions is, that the constrain about the concentrations does not allow two positive concentrations. Have a look at the following picture:

$\hskip3.5cm$enter image description here

The pH of a pure 0.0375 M sulfuric acid solution would be somewhat around 1.35. Changing the mixture by reducing the amount of sulfuric acid and increasing the amount of sulfurous acid leads to a rising pH. A pure 0.0375 M sulfurous acid solution would then have a pH of about 1.77.

But for concentrations below around $\mathrm{0.025~mol/L}$ the system is solvable, what can be seen in the next image. You could obtain the pH of 1.87 with a mixture of 0.00086859742 M $\ce{H2SO4}$ and 0.024131403 M $\ce{H2SO3}$.

$\hskip3.5cm$enter image description here

Now doing this with Octave

I never used Octave before, so there might be a better solution. But this is my code:

>> function y = Conc(x)
ka1 = 10^3; % wikipedia, sulfuric_acid
ka2 = 10^(-1.99); % wikipedia, sulfuric_acid
ka3 = 10^(-1.857); % wikipedia, sulfurous_acid
ka4 = 10^(-7.172); % wikipedia, sulfurous_acid
ctot = 0.025; % total sulfur concentration in mol/L
ph = 1.87; h = 10^-(ph);
oh = 10^(-14)/h;
xH2SO4 = (h * ka1 + 2 * ka1 * ka2)/(h^2 + h * ka1 + ka1 * ka2);
xH2SO3 = (h * ka3 + 2 * ka3 * ka4)/(h^2 + h * ka3 + ka3 * ka4);
y(1) = x(1) + x(2) - ctot;
y(2) = oh - h + x(1) * xH2SO4 + x(2) * xH2SO3;
>> [x, fval, info] = fsolve(@Conc, [1; 2])
x =
fval =
  6.3704e-011  2.8976e-011
info =  1

And as you can see, it gives exactle the same solutions as Mathematica. (fval is the value for y(1) and y(2) that are near zero, so everything is fine. info = 1 means, that everything the solution has converged.)

  • 1
    $\begingroup$ Yes, this is what I recognized as well. 1.87 is not accessible from the 0.0375 M system. Thank you for your excellent answer. $\endgroup$
    – JohnS
    Commented Aug 12, 2015 at 22:57

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.