6.10: Using Excel and R to Solve Equilibrium Problems
- Page ID
- 5720
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)In solving equilibrium problems we typically make one or more assumptions to simplify the algebra. These assumptions are important because they allow us to reduce the problem to an equation in x that we can solve by simply taking a square-root, a cube-root, or by using the quadratic equation. Without these assumptions, most equilibrium problems result in a cubic equation (or a higher-order equation) that is harder to solve. Both Excel and R are useful tools for solving such equations.
6.10.1 Excel
Excel offers a useful tool—the Solver function—for finding a polynomial equation’s chemically significant root. In addition, it is easy to solve a system of simultaneous equations by constructing a spreadsheet that allows you to test and evaluate multiple solutions. Let’s work through two examples.
Example 1: Solubility of Pb(IO3)2 in 0.10 M Pb(NO3)2
In our earlier treatment of this problem we arrived at the following cubic equation
\[4x^3+0.40x^2=2.5\times10^{-13}\]
where x is the equilibrium concentration of Pb2+. Although there are several approaches for solving cubic equations, none are computationally easy. One approach is to iterate in on the answer by finding two values of x, one of which leads to a result larger than 2.5×10–13 and one of which gives a result smaller than 2.5×10–13. Having established boundaries for the value of x, we then shift the upper limit and the lower limit until the precision of our answer is satisfactory. Without going into details, this is how Excel’s Solver function works.
To solve this problem, we first rewrite the cubic equation so that its right-side equals zero.
\[4x^3+0.40x^2-2.5\times10^{-13}=0\]
Next, we set up the spreadsheet shown in Figure 6.17a, placing the formula for the cubic equation in cell B2, and entering our initial guess for x in cell B1. We expect x to be small—because Pb(IO3)2 is not very soluble—so setting our initial guess to 0 seems reasonable. Finally, we access the Solver function by selecting Solver... from the Tools menu, which opens the Solver Parameters window.
Figure 6.17 Spreadsheet demonstrating the use of Excel’s Solver function to find the root of a cubic equation. The spreadsheet in (a) shows the cubic equation in cell B2 and the initial guess for the value of x in cell B1; Excel replaces the formula with its equivalent value. The spreadsheet in (b) shows the results of running Excel’s Solver function.
To define the problem, place the cursor in the box for Set Target Cell and then click on cell B2. Select the Value of: radio button and enter 0 in the box. Place the cursor in the box for By Changing Cells: and click on cell B1. Together, these actions instruct the Solver function to change the value of x, which is in cell B1, until the cubic equation in cell B2 equals zero.
Before we actually solve the function, we need to consider whether there are any limitations for an acceptable result. For example, we know that x cannot be smaller than 0 because a negative concentration is not possible. We also want to ensure that the solution’s precision is acceptable. Click on the button labeled Options... to open the Solver Options window. Checking the option for Assume Non-Negative forces the Solver to maintain a positive value for the contents of cell B1, meeting one of our criteria. Setting the precision takes a bit more thought. The Solver function uses the precision to decide when to stop its search, doing so when
\[|expected\;value-calculated\;value|\times100\leq precision\,(\%)\]
where expected value is the target cell’s desired value (0 in this case), calculated value is the function’s current value (cell B1 in this case), and precision is the value we enter in the box for Precision. Because our initial guess of x = 0 gives a calculated result of 2.5×10–13, accepting the Solver’s default precision of 1×10–6 stops the search after one cycle. To be safe, set the precision to 1×10–18. Click OK and then Solve. When the Solver function finds a solution, the results appear in your spreadsheet (see Figure 6.17b). Click OK to keep the result, or Cancel to return to the original values. Note that the answer here agrees with our earlier result of 7.91×10–7 M for the solubility of Pb(IO3)2.
Note
Be sure to evaluate the reasonableness of Solver’s answer. If necessary, repeat the process using a smaller value for the precision.
Example 2: pH of 1.0 M HF
In developing our solution to this problem we began by identifying four unknowns and writing out the following four equations.
\[K_\mathrm a=\mathrm{\dfrac{[H_3O^+][F^-]}{[HF]}}=6.8\times10^{-4}\]
\[K_\mathrm w=\mathrm{[H_3O^+][OH^-]}=1.00\times10^{-14}\]
\[C_\mathrm{HF}=\mathrm{[HF]+[F^-]}\]
\[\mathrm{[H_3O^+]=[OH^-]+[F^-]}\]
From this point, we made two assumptions, simplifying the problem to one that was easy to solve.
\[\mathrm{[H_3O^+]}=\sqrt{K_\mathrm aC_\mathrm{HF}}=\sqrt{(6.8\times10^{-4})(1.0)}=2.6\times10^{-2}\]
Although we did not note this at the time, without making assumptions the solution to our problem is a cubic equation
\[\mathrm{[H_3O^+]^3}+K_\mathrm a\mathrm{[H_3O^+]^2}-(K_\mathrm aC_\mathrm{HA}+K_\mathrm w)\mathrm{[H_3O^+]}-K_\mathrm aK_\mathrm w=0\tag{6.64}\]
that we can solve using Excel’s Solver function. Of course, this assumes that we successfully complete the derivation!
Another option is to use Excel to solve the equations simultaneously by iterating in on values for [HF], [F−], [H3O+], and [OH−]. Figure 6.18a shows a spreadsheet for this purpose. The cells in the first row contain our initial guesses for the equilibrium pH. Using the ladder diagram in Figure 6.14, choosing pH values between 1 and 3 seems reasonable. You can add additional columns if you wish to include more pH values. The formulas in rows 2–5 use the definition of pH to calculate [H3O+], Kw to calculate [OH−], the charge balance equation to calculate [F−], and Ka to calculate [HF]. To evaluate the initial guesses, we use the mass balance expression for HF, rewriting it as
\[\mathrm{[HF]+[F^-]-\mathit C_{HF}=[HF]+[F^-]-1.0=0}\]
and entering it in the last row. This cell gives the calculation’s error.
Figure 6.18b shows the actual values for the spreadsheet in Figure 6.18a.The negative value in cells B6 and C6 means that the combined concentrations of HF and F– are too small, and the positive value in cell D6 means that the combined concentrations are too large. The actual pH, therefore, must lie between 2.00 and 1.00. Using these pH values as new limits for the spreadsheet’s first row, we continue to narrow the range for the actual pH. Figure 6.18c shows a final set of guesses, with the actual pH falling between 1.59 and 1.58. Because the error for 1.59 is smaller than that for 1.58, we will accept a pH of 1.59 as the answer. Note that this is an agreement with our earlier result.
Figure 6.18 Spreadsheet demonstrating the use of Excel to solve a set of simultaneous equations. The spreadsheet in (a) shows the initial guess for [H3O+] in the first row, and the formulas that you must enter in rows 2–6. Enter the formulas in cells B2–B6 and then copy and paste them into the appropriate cells in the remaining columns. As shown in (b), Excel replaces the formulas with their equivalent values. The spreadsheet in (c) shows the results after our final iteration. See the text for further details.
Note
You also can solve this set of simultaneous equations using Excel’s Solver function. To do so, create the spreadsheet in Figure 6.18a, but omit all columns other than A and B. Select Solver... from the Tools menu and define the problem by using B6 for Set Target Cell and setting its desired value to 0, and selecting B1 for By Changing Cells:. You may need to play with the Solver’s options to find a suitable solution to the problem, and it is wise to try several different initial guesses.
The Solver function works well for relatively simple problems, such as finding the pH of 1.0 M HF. As problems become more complex—such as solving an equilibrium problem with lots of unknowns—the Solver function becomes less reliable in finding a solution.
Practice Exercise 6.14
Using Excel, calculate the solubility of AgI in 0.10 M NH3 without making any assumptions. See our earlier treatment of this problem for the relevant equilibrium reactions and constants.
Click here to review your answer to this exercise.
6.10.2 R
R has a simple command—uniroot—for finding the chemically significant root of a polynomial equation. In addition, it is easy to write a function to solve a set of simultaneous equations by iterating in on a solution. Let’s work through two examples.
Example 1: Solubility of Pb(IO3)2 in 0.10 M Pb(NO3)2
In our earlier treatment of this problem we arrived at the following cubic equation
\[4x^3+0.40x^2=2.5\times10^{-13}\]
where x is the equilibrium concentration of Pb2+. Although there are several approaches for solving cubic equations, none are computationally easy. One approach to solving the problem is to iterate in on the answer by finding two values of x, one of which leads to a result larger than 2.5×10–13 and one of which gives a result smaller than 2.5×10–13. Having established boundaries for the value of x, we then shift the upper limit and the lower limit until the precision of our answer is satisfactory. Without going into details, this is how the uniroot command works.
The general form of the uniroot command is
uniroot(function, lower, upper, tol)
where function is an object containing the equation whose root we are seeking, lower and upper are boundaries for the root, and tol is the desired accuracy for the root. To create an object containing the equation, we must rewrite it so that its right-side equals zero.
\[4x^3+0.40x^2-2.5\times10^{-13}=0\]
Next, we enter the following code, which defines our cubic equation as a function with the name eqn.
> eqn = function(x) 4*x^3 + 0.4*x^2 – 2.5e–13
Because our equation is a function, the uniroot command can send a value of x to eqn and receive back the equation’s corresponding value.
Note
For example, entering
> eqn(2)
passes the value x = 2 to the function and returns an answer of 33.6.
Finally, we use the uniroot command to find the root.
> uniroot(eqn, lower = 0, upper = 0.1, tol = 1e–18)
We expect x to be small—because Pb(IO3)2 is not very soluble—so setting the lower limit to 0 is a reasonable choice. The choice for the upper limit is less critical. To ensure that the solution has sufficient precision, the tolerance should be smaller than the expected root. Figure 6.19 shows the resulting output. The value $root is the equation’s root, which is in good agreement with our earlier result of 7.91×10–7 for the molar solubility of Pb(IO3)2. The other results are the equation’s value for the root, the number of iterations needed to find the root, and the root’s estimated precision.
Figure 6.19 The summary of R’s output from the uniroot command. See the text for a discussion of how to interpret the results.
Example 2: pH of 1.0 M HF
In developing our solution to this problem we began by identifying four unknowns and writing out the following four equations.
\[K_\mathrm a=\mathrm{\dfrac{[H_3O^+][F^-]}{[HF]}}=6.8\times10^{-4}\]
\[K_\mathrm w=\mathrm{[H_3O^+][OH^-]}=1.00\times10^{-14}\]
\[C_\mathrm{HF}=\mathrm{[HF]+[F^-]}\]
\[\mathrm{[H_3O^+]=[OH^-]+[F^-]}\]
From this point, we made two assumptions, simplifying the problem to one that was easy to solve.
\[\mathrm{[H_3O^+]}=\sqrt{K_\mathrm aC_\mathrm{HF}}=\sqrt{(6.8\times10^{-4})(1.0)}=2.6\times10^{-2}\]
Although we did not note this at the time, without making assumptions the solution to our problem is a cubic equation
\[\mathrm{[H_3O^+]^3+\mathit K_a[H_3O^+]^2-(\mathit K_aC_{HA}+\mathit K_w)[H_3O^+]-\mathit K_a\mathit K_w=0}\]
that we can solve using the uniroot command. Of course, this assumes that we successfully complete the derivation!
Another option is to use write a function to solve simultaneously the four equations for the variables [HF], [F−], [H3O+], and [OH−]. Here is the code for this function, which we will call eval.
> eval = function(pH){
+ h3o =10^–pH
+ oh = 1e–14/h3o
+ hf = (h3o*f)/6.8e–4
+ error = hf + f – 1
+ output = data.frame(pH, error)
+ print(output)
+ }
Note
The open { tells R that we intend to enter our function over several lines. When we press enter at the end of a line, R changes its prompt from > to +, indicating that we are continuing to enter the same command. The close } on the last line indicates that we are done entering the function.
The command data.frame combines two or more objects into a table.
You can adapt this function to other problems by changing the variable you pass to the function and the equations you include within the function.
Let’s examine more closely how this function works. The function accepts a guess for the pH and uses the definition of pH to calculate [H3O+], Kw to calculate [OH−], the charge balance equation to calculate [F−], and Ka to calculate [HF]. The function then evaluates the solution using the mass balance expression for HF, rewriting it as
\[\mathrm{[HF]+[F^-]-\mathit C_{HF}=[HF]+[F^-]-1.0=0}\]
The function then gathers together the initial guess for the pH and the error and prints them as a table.
The beauty of this function is that the object we pass to it, pH, can contain many values, allowing us to easily search for a solution. Because HF is an acid, we know that the solution must be acidic. This sets a lower limit of 7 for the pH. We also know that the pH of 1.0 M HF can be no larger than 1.0 M as this would be the upper limit if HF was a strong acid. For our first pass, let’s enter the following code
> pH = c(7, 6, 5, 4, 3, 2, 1)
> func(pH)
which varies the pH within these limits. The result, which is shown in Figure 6.20a, indicates that the pH must be less than 2 and greater than 1 because it is in this interval that the error changes sign.
For our second pass, we explore pH values between 2.0 and 1.0 to further narrow down the problem’s solution.
> pH = c(2.0, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0)
> func(pH)
The result in Figure 6.20b show that the pH must be less than 1.6 and greater than 1.5. A third pass between these limits gives the result shown in Figure 6.20c, which is consistent with our earlier result of a pH 1.59.
Figure 6.20 The output of three iterations to find the pH for a solution of 1.0 M HF. The results are for pH values between (a) 7 and 0, (b) 2.0 and 1.0, and (c) 1.60 M and 1.50. The columns labeled “error” show an evaluation of the mass balance equation for HF, with positive values indicating that the pH is too low and negative values indicating that the pH is too high.
Practice Exercise 6.15
Using R, calculate the solubility of AgI in 0.10 M NH3 without making any assumptions. See our earlier treatment of this problem for the relevant equilibrium reactions and constants.
Click here to review your answer to this exercise.