We model the titration of a weak triprotic acid by a strong base using an arc-length continuation method in Mathematica. The acid is characterized by three dissociation constants corresponding to pKa1=3, pKa2=6, and pKa3=9 (see Ref. 1). During the titration, the solution undergoes successive deprotonation reactions, producing the species
H
3
A,
H
2
-
A
, H
2-
A
, and
3-
A
. ​A very common real-world example of a triprotic acid is Phosphoric acid. Indeed, Phosphoric acid has three acidic protons that dissociate step by step:
H
3
PO4 ⇌
+
H
+
H
2
PO
-
4
​
H
2
PO
-
4
⇌
+
H
+
H
PO
2-
4
​
H
PO
2-
4
⇌
+
H
+ PO
3-
4
​Its approximate dissociation constants are: pKa1≈2.15, pKa2≈7.20, pKa3≈12.35. This acid is widely used in: fertilizers, food and beverages (especially soft drinks), rust removal, water treatment, biochemical buffer solutions.
The mathematical model is based on chemical equilibrium relations, mass conservation, and electroneutrality. The nonlinear algebraic system describing the titration curve may become difficult to solve near equivalence regions because of sharp pH variations and multiple buffering zones. To improve numerical robustness, the arc-length continuation method is employed. Instead of solving independently for each added-base fraction, the continuation algorithm follows the equilibrium solution branch progressively using predictor-corrector steps. This approach ensures stable convergence even near steep transitions and equivalence points.
The simulation provides both the pH titration curve and the distribution diagram of the chemical species as a function of the normalized titrant amount x defined as moles of base added over initial moles of acid.
​
The results clearly show the three buffer regions associated with the three pKa​ values, as well as the dominance regions of each species during titration.
The method demonstrates the effectiveness of arc-length continuation for tracking highly nonlinear equilibrium systems encountered in acid-base chemistry.
Acid concentration mol/l
Ca=0.01;
NaOH concentration mol/l
Cb=0.1;
Initial acid volume l
Va=0.05;
Dissociation constants
pKa1=3;​​pKa2=6;​​pKa3=9;​​​​Ka1=10^-pKa1;​​Ka2=10^-pKa2;​​Ka3=10^-pKa3;
Ionic product of water (also called the water dissociation constant)
Kw=10^-14;
Initial number of moles of acid present before titration begins
na0=CaVa;​​
Species distribution functions
alpha0[h_]:=h^3/(h^3+Ka1h^2+Ka1Ka2h+Ka1Ka2Ka3);​​alpha1[h_]:=Ka1h^2/(h^3+Ka1h^2+Ka1Ka2h+Ka1Ka2Ka3);​​alpha2[h_]:=Ka1Ka2h/(h^3+Ka1h^2+Ka1Ka2h+Ka1Ka2Ka3);​​alpha3[h_]:=Ka1Ka2Ka3/(h^3+Ka1h^2+Ka1Ka2h+Ka1Ka2Ka3);
Total volume
Vtot[x_]:=Va+(xna0)/Cb;
charge balance equation
x=nb/na0; x is defined as moles of base added over initial moles of acid.
f[h_,x_]:=Module[{Ct,na,nb},na=na0;​​nb=xna0;​​Ct=na/Vtot[x];​​h+(nb/Vtot[x])-Kw/h-Ct(alpha1[h]+2alpha2[h]+3alpha3[h])];
Arc-length continuation
Unknowns={h,x}
ds=0.01;(*arcstep*)​​nSteps=400;
(*Initialpoint*)​​x0=0.0;​​h0=10^-2.6;​​​​sol={{x0,h0}};​​​​x=x0;​​h=h0;
Continuation loop; three steps :
(1) computes a tangent vector
(2) performs a predictor step along the tangent, then
(3) applies a corrector step using Newton iterations
Do[​​(*Jacobiancomponents*)​​fx=D[f[hvar,xvar],xvar]/.{hvar->h,xvar->x};​​fh=D[f[hvar,xvar],hvar]/.{hvar->h,xvar->x};​​​​(*Tangentvector*)​​tx=1;​​th=-fx/fh;​​​​(*Normalize*)​​norm=Sqrt[tx^2+th^2];​​tx=tx/norm;​​th=th/norm;​​​​(*Predictor*)​​xPred=x+dstx;​​hPred=h+dsth;​​​​(*Correctorsystem*)solNew=FindRoot[{f[hnew,xnew]==0,(xnew-xPred)tx+(hnew-hPred)th==0},{xnew,xPred},{hnew,hPred},MaxIterations->100];​​​​x=xnew/.solNew;​​h=hnew/.solNew;​​​​AppendTo[sol,{x,h}],{nSteps}]
Extract data
xVals=sol[[All,1]];​​hVals=sol[[All,2]];
pH
Distribution fractions
Distribution diagram

References

[1] https://ressources.univ-lemans.fr/AccesLibre/UM/Pedago/chimie/01/05-Chimie-en-solution/res/Acidimetrie-acides-dotriac.pdf

CITE THIS NOTEBOOK

Titration of a triprotic acid by strong base​
by Housam Binous​
Wolfram Community, STAFF PICKS, May 25, 2026
https://community.wolfram.com/groups/-/m/t/3721801