suppress that Solve: message

Wrap it your Solve[...] expression with Quiet. So Quiet[Solve[...]] or Solve[...]//Quiet

why, after making a change to the code, I have to Evaluate the notebook twice for the changes to be reflected?

I'm just guessing, but it probably has to do with all these symbols being global and therefore the sequence of computation is important.

So, for example if at the top of my notebook a evaluate a=7 and then later evaluate b=a, at that point we have "facts" that a is 7 and b is 7. Now if you evaluate a=13, you've only changed one of those facts, i.e. a is now 13 but b is still 7. If instead you had written b := a, then b would "pick up" the value of a each time it is evaluated.

But I haven't gone through your notebook recently, so that's just a guess.

It would be nice if I could somehow suppress that Solve: message

Use exact numerical values

x = Solve[Rationalize[makeRelationEquation[Ka[pKa1], Kb[pKb], Ca1, Cb1, Va], 0], H];

Use Quiet with caution. It will suppress all messages.

With exact input the result is also exact. In this case Root objects. Instead of the complicated nested Part you can do this

sol = H /. x

Which will return a list of the 4 solutions in sol. To get the numerical (inexact values)

sol // N

To get exact representations of the Root objects

sol // ToRadicals

Paul,

Please post complete working code. relationEquation2 refers to symbols Ca and Cb that are not defined anywhere earlier in this thread. Here is what I evaluated

Kw = 1.00*10^-14;
Ka[pKa_] := 10^-pKa;
Kb[pKb_] := 10^-pKb;
pKa1 := 4.757;
pKb := 0.63;
Ca1 := 0.1;
Cb1 := 0.1;
Va := 10;
Vb := 10;
relationEquation = 
  Vb == Va*((Ca/(1 + H/Ka)) - H + Kw/H)/((Cb/(1 + Kw/(H*Kb))) + H - Kw/H);
makeRelationEquation[Ka_, Kb_, Ca_, Cb_, Va_] = relationEquation;
x = Solve[Rationalize[makeRelationEquation[Ka[pKa1], Kb[pKb], Ca1, Cb1, Va], 0], H];

Extract the solutions, select the positive one and convert to a numerical value

H /. x // Select[Positive] // First // N
(* 2.06037*10^-9 *)

Exact value

H /. x // Select[Positive] // First // ToRadicals
(* Long expression *)

Another option is to use SolveValues, then ReplaceAll is not needed.

x = SolveValues[
  Rationalize[makeRelationEquation[Ka[pKa1], Kb[pKb], Ca1, Cb1, Va], 0], H]
x // Select[Positive] // First // N

/. is ReplaceAll.

I believe I may have a solution. From a chemistry perspective, the only valid solutions are those where H > 0. I believe that only one solution will meet the criterion H > 0 and although I could be wrong I propose to use that as a working hypothesis. In which case I can do this

x = Solve[makeRelationEquation[Ka[pKa1],Kb[pKb],Ca1,Cb1,Va],H,Assumptions->H>0];
x
pH = -Log10[Part[Part[Part[x,1],1],2]]

Do you have any other ideas? I have tested this with these values:

Kw = 1.00*10^-14;
Ka[pKa_] := 10^-pKa;
Kb[pKb_] := 10^-pKb;
pKb := 0.2;
Ca1 := 0.1;
Cb1 := 0.1;
Va := 10;
Vb := 10;
pKa := 0.63 and pKa := 4.757

Hi Eric, I am trying to adapt some of the code that you have previously provided to give me H (and hence pH) for just one value of Vb and this is what I have so far:

Kw = 1.00*10^-14;
Ka[pKa_] := 10^-pKa;
Kb[pKb_] := 10^-pKb;
pKa1 := 4.757;
pKb := 0.63;
Ca1 := 0.1;
Cb1 := 0.1;
Va := 10;
Vb := 10;
relationEquation=Vb==Va*((Ca/(1+H/Ka))-H+Kw/H)/((Cb/(1+Kw/(H*Kb)))+H-Kw/H);
makeRelationEquation[Ka_,Kb_,Ca_,Cb_,Va_]=relationEquation;
x = Solve[makeRelationEquation[Ka[pKa1],Kb[pKb],Ca1,Cb1,Va],H];
x
pH = -Log10[2.06037 * 10^-9] 

This gives the following:

Obviously I had to select which of the four values in x was the sensible one and I wondered if there is a better, more direct way of getting the result. Thanks.

I don't know if you are still keeping an eye on this discussion, Eric, but if you are I thought you might be interested in what I have finally come up with. If you, or anybody else, is interested I will happily attach my notebook (such as it is).

It's a first attempt and I am sure that you could pick it to pieces ..anyway you can find it here Trying to add a legend when using CombinePlots I couldn't have got even this far without your considerable help, so thank you once again.

If this is all true, - it is the truth, the whole truth and nothing but the truth

Yes, exactly, we can let the user manipulate Ka, Kb, Ca, Cb and Va

we want to show this in either orientation (i.e. for a target pH we want to see what Vb volume we need, or alternatively for a given Vb volume we want to know the pH we achieve

Yes again, but in both cases I would prefer pH on the y-axis

I had already sorted default values in Manipulate and that (and other things) are on the last notebook I attached and I did follow up on all the things you sent and referred to.

Sorry for the confusion - "we don't need to know anything about Ka or Kb" applies in the laboratory situation when we are dealing with real solutions and glassware. Again, we are only solving for Ca in the laboratory.

So to clarify (and this will make more sense when you look at my latest version of the notebook which I have attached and to which I have added a number of enhancements)

Kw is always a constant 1.0 E-14

Ka, Kb are dependent on the acid and base respectively and will always be known

Ca, Cb are the concentrations of the acid and base respectively - and in the case of our Mathematica approach will be determined by the user

Va is the volume of acid to which we are adding the base - the user sets this value, usually to 5, 10, 20 or 25

What we are trying to achieve is a plot of pH vs Vb. We can do this in one of two ways;

The first way (which we have been working on) is to calculate Vb from pH. The resulting graph has Vb on the y-axis and pH on the x-axis - I don't know if the axes can be transposed, but it would be nice if they could

The second way (which is what we want to do next) is to calculate pH from Vb and then plot pH (y-axis) against Vb (x-axis). You said that you might find a way " I think we might be able to coerce Mathematica into solving the big hairy formula for pH." It might be a good idea to do this in a separate notebook.

I hope all this now makes sense but, if not, please let me know.

Attachments:

Titration (Rev 2).nb

Okay, I think I understand that. What I'm confused about is what you're trying to find. Previously, it sounded like you were solving for pH, but now it sounds like you're solving for Ca. And when we started, we were solving for Vb. I guess I'm also confused about this: "we don't need to know anything about Ka or Kb". Either we have values for Ka & Kb, or we don't. They seem to be present in the formula you referenced, so I don't know how to set this up for the case where we don't know those values.

Here's what it sounds like to me. We have some acid solution. We must know some of its properties, but I'm confused about which ones. What we don't know but would like to know is its concentration, Ca. We have some base solutions, and we know some of their properties. If we can get to a neutral solution (I assume that's the "endpoint", but maybe it's some other pH) by adding a known amount of the base, we can calculate Ca with this complicated formula. So, we're solving for Ca, and we need a formula where Ca is defined in terms of all of these other independent variables. It sounds like the goal is to solve that complicated formula for Ca. Am I on the right track?

Sure. Sorry about the units (can't do superscripts here) Ca and Cb are the concentrations in mol dm-3 of an acid and a base Va and Vb are the volumes in cm3 of the acid and base Ka and Kb are the dissociation constants of the acid and base in mol dm-3 Kw is the ionic product of water in mol2 dm-6

The equation comes from an article in the journal School Science Review

Greatorex, D., 'Exact algorithms for acid/base titration curves', S.S.R. 1979, 214, 61, 94

I eventually got hold of a copy (attached) with help of a very kind ladyn at the Association for Science Education (the publishers of S.S.R.)

Attachments:

Greatorex, D., '...pdf

Vis-a-vis the divide by zero range error, it's a shame that instead of (for example) {Kb,0,1} which means 0 <= Kb <= 1 there could be way of stating the range as 0 < Kb <= 1

2.

...what I really want to do is calculate pH for a range of values of Vb...

If you can rewrite the equation in terms of ph = <some complicated expression>, then sure, that's doable. I'm not going to try to do it myself :). Then we'll just write a definition like ph[<arguments go here>]:=<new definition goes here>.

3B.

...Also can the number of decimal places in the result be specified...

Look at the various *Form methods: http://reference.wolfram.com/language/guide/DisplayOfNumbers.html

Try NumberForm, for example.

5.

is it possible to have the table output arranged as two columns

TableForm[Table[{x, Vb[x]}, {x, 3, 12.5, 0.2}]]

If you're going to want to do further processing on that data, then I suggest you first save the data in a variable:

dataPairs = Table[{x, Vb[x]}, {x, 3, 12.5, 0.2}];
TableForm[dataPairs]

The various *Form functions are just special wrappers that the front end knows how to interpret to give certain visualizations. You rarely want to use form-wrapped expressions in computations. In other words, you probably do NOT want to do this:

dataPairs = TableForm[Table[{x, Vb[x]}, {x, 3, 12.5, 0.2}]]

TableForm has a bunch of options if you want to tweak further.

The underscores are a pattern matching construct. Also, the := is SetDelayed, which is a common way to define a function.

http://reference.wolfram.com/language/ref/Blank.html http://reference.wolfram.com/language/ref/SetDelayed.html https://www.wolfram.com/language/elementary-introduction/2nd-ed/39-immediate-and-delayed-values.html https://www.wolfram.com/language/fast-introduction-for-programmers/en/function-definitions/