March 28, 2024, 03:38:30 PM
Forum Rules: Read This Before Posting


Topic: Titration curve (MATLAB simulation)  (Read 10932 times)

0 Members and 1 Guest are viewing this topic.

Offline PoetryInMotion

  • Regular Member
  • ***
  • Posts: 49
  • Mole Snacks: +1/-0
Titration curve (MATLAB simulation)
« on: June 29, 2015, 11:41:56 AM »
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 C0 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.

Code: [Select]
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.
« Last Edit: June 29, 2015, 12:04:02 PM by PoetryInMotion »
Undergraduate student majoring in chemistry and mathematics. Former IChO participant.

Offline mjc123

  • Chemist
  • Sr. Member
  • *
  • Posts: 2048
  • Mole Snacks: +296/-12
Re: Titration curve (MATLAB simulation)
« Reply #1 on: June 29, 2015, 11:59:58 AM »
Not sure how your code works, but when you write
K = exp(-4.76)
does that mean 10-4.76 or e-4.76?
If the latter, then your pKa, and pH at half titre, would be 2.1 - which looks close to what the graph says.

Offline Arkcon

  • Retired Staff
  • Sr. Member
  • *
  • Posts: 7367
  • Mole Snacks: +533/-147
Re: Titration curve (MATLAB simulation)
« Reply #2 on: June 29, 2015, 12:10:52 PM »
I don't understand what you claim to be missing, looking at your examples, there isn't much of a graphical structure associated with half-titration point.
Hey, I'm not judging.  I just like to shoot straight.  I'm a man of science.

Offline PoetryInMotion

  • Regular Member
  • ***
  • Posts: 49
  • Mole Snacks: +1/-0
Re: Titration curve (MATLAB simulation)
« Reply #3 on: June 29, 2015, 12:19:03 PM »
Not sure how your code works, but when you write
K = exp(-4.76)
does that mean 10-4.76 or e-4.76?
If the latter, then your pKa, and pH at half titre, would be 2.1 - which looks close to what the graph says.
Gah! Such an embarrassing mistake :D

Lesson of today: never underestimate your own stupidity when debugging your code.

New result:



I am incredibly thankful.
Undergraduate student majoring in chemistry and mathematics. Former IChO participant.

Offline samcarano

  • Very New Member
  • *
  • Posts: 2
  • Mole Snacks: +0/-0
Re: Titration curve (MATLAB simulation)
« Reply #4 on: April 03, 2019, 05:00:01 PM »
I am trying to work through the derivations but I don't understand where the last two of the 4 equations come from.
Would somebody mind explaining to me where Ca0 and Cb0 are obtained from? 

The equations I see besides the ones mentioned here that seem relevant:
[HA]-->[H+]+[A-]
[HA]+[OH]-->[A-]+[H+]

I would think Ca0 the same thing as [HA] and Cb0 the same as [OH]

Clearly that is not true here but I can't seem to figure what I am not getting.

Any help would be much apprecitaed,

Thanks
« Last Edit: April 03, 2019, 05:26:18 PM by samcarano »

Online Borek

  • Mr. pH
  • Administrator
  • Deity Member
  • *
  • Posts: 27633
  • Mole Snacks: +1799/-410
  • Gender: Male
  • I am known to be occasionally wrong.
    • Chembuddy
Re: Titration curve (MATLAB simulation)
« Reply #5 on: April 03, 2019, 05:56:27 PM »
Would somebody mind explaining to me where Ca0 and Cb0 are obtained from?

They aren't "obtained", they are defined.

They are analytical concentrations. When you put acid HA in the solution it dissociates, you have HA and A- present - and even if you don't know the dissociation extent it is obvious that the sum of both forms must be identical to the amount of acid put in the solution. Say, you put 1 mole of acid in 1000 mL - no matter what the Ka is, [HA]+[A-]=1.
ChemBuddy chemical calculators - stoichiometry, pH, concentration, buffer preparation, titrations.info

Offline samcarano

  • Very New Member
  • *
  • Posts: 2
  • Mole Snacks: +0/-0
Re: Titration curve (MATLAB simulation)
« Reply #6 on: April 03, 2019, 06:25:44 PM »
That makes some sense.  Thank you.  So you are adding the acid concentration and the acid conjugate concentration.

Cb0 still confuses me:

Why is Cb0= [OH-]-[H+]+[A-]?
I would think it would be just Cb0= [OH-]+[H+] to represent the total base in the solution.  Where does the [A-] come from and why not add the base and its conjugate again?

I might be getting throw off by what they poster means by "imaginary concentration"

Thanks for the help

Online Borek

  • Mr. pH
  • Administrator
  • Deity Member
  • *
  • Posts: 27633
  • Mole Snacks: +1799/-410
  • Gender: Male
  • I am known to be occasionally wrong.
    • Chembuddy
Re: Titration curve (MATLAB simulation)
« Reply #7 on: April 04, 2019, 03:01:37 AM »
Yeah, it can be confusing. Cb0 equation combines three things - mass balance for base, assumption that base is strong and fully dissociated, and charge balance.

BOH  :rarrow: B+ + OH-

Start with

Cb0=[B+] + [BOH]

write charge balance and see what you get from these.
ChemBuddy chemical calculators - stoichiometry, pH, concentration, buffer preparation, titrations.info

Sponsored Links