I have written a small MATLAB script that simulates a titration curve for the titration of a weak acid with NaOH (or similar) as titrator.

The result when I simulate the titration of 0.5 M acetic acid (pKa=4.76) with 0.5 M NaOH is this:

What bothers me is that I do not see any clear half-titration point. There ought to be one, according to with what seem to be experimental titration curves

here and

here. Hopefully someone here can help me understand why my script does not work. I have tried my best to eliminate approximation (for instance that the volume of the solution is approximately constant or that the self-dissociation of water does not affect the pH) but perhaps physical chemistry above my capacity is needed to get enough accuracy to see the half-titration point.

Below I include the code and the mathematical background, so that you get a better picture of what I have done this far.

**Mathematical principle.**What my script does, is that it after each small addition of base, compute the pH in the solution. Both the protolysis equilibrium and the self-dissociation of water is taken into account. More precisely, I solve the following system numerically:

[tex]\begin{cases}K[HA]=[H^+][A^-]\\ Kw=[H^+][OH^-]\\ C_{a0} = [HA] + [A-]\\ C_{b0} = [OH^-]-[H^+] + [A-]\,,\end{cases}[/tex]

where K is the acid constant for my acid, Kw is the self-dissociation constant for water and C

_{0} represents the (imaginary) concentration of acid and base before the system is allowed to equilibriate.

My best guess for why

Substitutions give rise to

[tex]C_{a0}=\frac{[H^+][A^-]}{[K]}+[A^-]=(\frac{[H^]}{K}+1)(C_{b0}+[H^+]-\frac{K_w}{[H^+]})\,,[/tex]

[tex]K[H^+]C_{a0}=([H^+]+K)(C_{b0}[H^+]+[H^+]^2-K_w)\,,[/tex]

and finally

[tex][H^+]^3 + (K+C_{b0})[H^+]^2 + (C_{b0}K - C_{a0}K - K_w)[H^+] -KK_w=0\,,[/tex]

which is a polynomial which can easily be solve with a built in function in MATLAB.

**My code.**`Kw = 1e-14;`

% Initial concentration of acid

C_acid = 0.5; % mol/dm^3

% Molarity of added base

C_base = 0.5; % mol/dm^3

% Initial volume of acid

vol = 25.0; % ml

% Stop graph when x ml base is added

max_added = 50.0; % ml

% Acid constant

K = exp(-4.76);

% Graph settings

dx = 0.0001;

x = 0.0:dx:max_added;

% pH = zeros(10 * max_added)

% Calculate the initial number of moles of acid

n_a = C_acid * vol * 0.001;

for i = 1:1/dx*max_added+1

% Calculate the number of moles of base added

n_b = x(i) * C_base * 0.001;

% The total volume

vol_tot = vol + x(i);

% Initial concentrations of acid and base

C_a0 = n_a / vol_tot * 1000;

C_b0 = n_b / vol_tot * 1000;

% Compute the concentration of hydrogen ions

% Solve the polynomial equation

R = roots([1, (K+C_b0), (C_b0*K - C_a0*K - 1e-14), -K*1e-14]);

% Remove all chemically irrelevant solutions

R = R(R==real(R) & R>0 & R<C_a0);

% If there are more than one chemically relevant solution, display the number

if not(length(R) == 1)

disp(length(R))

end

H = min(R(R==real(R) & R>0 & R<C_a0));

% Calculate the pH value

pH(i) = -log10(H);

end

% Plot it

plot(x,pH,'red')

xlabel('TIllsatt volym bas (ml)')

%set(gca,'XTick',0:1:max_added);

%set(gca,'YTick',round(pH(1)):1:round(pH(length(pH)))+1);

%set(gca,{'xminorgrid' 'yminorgrid'}, {'on' 'on'});

ylabel('pH')

title('Titreing av ättiksprit')

Edit: For some reason, the LaTeX tool does not seem to work for me.