kill(all)$

Maxima : intro computation ,base and acid, 'part :17b,version 1.0'
peter.vlasschaert@gmail.com, 10/01/2019

load packages
============

load(lrats)$
load(facexp)$

definition : log10(x)
================

log10(x):= log(x)/log(10);

How to find inflection points from titration curve.
--------------------------------------------------------------------------

titration of triprotic acid with strong base (NaOH).
=========================================

-----------------------------------------------------------------------------------------------
given : triprotic acid (H3PO4)  conc = C_Z  = 0.1 M,(vessel), V_Z
              strong   base (NaOH  ) conc = C_B  = 0.1 M,(buret ),  V_B
-----------------------------------------------------------------------------------------------
V_Z  : volume 'vessel'
V_B : volume 'buret  '
during titration : volume  added = V_B + V_Z
-----------------------------------------------------------------------------------------------
strong base : NaOH  -> Na(+) + OH(-) ' full dissociation '
-----------------------------------------------------------------------------------------------

H3A  , tri   , acid   => triprotic acid     (n=3)
==================================
A = PO4

distribution function : (n=3)
=======================
H3PO4 :  α(i)  , i = 0,1,2,...,n
-----------------------------------------
K_a1 : 10^(-2.23)
K_a2 : 10^(-7.21)
K_a3 :10^(-12.32)
K_w : 10^(-14.00)
-----------------------------------------

rem : [H3A] : [H3A], [H2A(-)]:[H2A],[HA(2-)]:[HA],[A(3-)]:[A]

p1:C_Z;

p1a:α(0)= H3A/p1;
p1b:α(1)= H2A/p1;
p1c:α(2)= HA/p1;
p1d:α(3)= A/p1;

p2:H^3+K_a1*H^2+K_a1*K_a2*H+K_a1*K_a2*K_a3;

p2a:α(0)= H^3/p2;
p2b:α(1)= K_a1*H^2/p2;
p2c:α(2)= K_a1*K_a2*H/p2;
p2d:α(3)=K_a1*K_a2*K_a3/p2;

concentration : C_Z
============

p3:C_Z = H3A+H2A+HA+A;

K_w : dissociation constant   'water'
=====

p3a: OH = K_w/H;

species (reaction): with charge  =>  charge balance:
===========================================

-----------------------------------------------------------------------------
[H(+)] +[Na(+)] = [H2A(-)] +[HA(2-)] +[A(3-)] + [OH(-)]
=> see 'above'
-----------------------------------------------------------------------------

p4:H+Na = H2A+2*HA+3*A+OH;

combination: α(i) and change volume during titration
============================================

1e) without change of volume
--------------------------------------------

p4a:solve(p1b,H2A)[1];
p4b:solve(p1c,HA)[1];
p4c:solve(p1d,A)[1];

2e) with change of volume
-----------------------------------------

p5:C_Z = C_Z*V_Z/(V_Z+V_B);

p5a:subst([p5],p4a);
p5b:subst([p5],p4b);
p5c:subst([p5],p4c);

p5d:Na=C_B*V_B/(V_Z+V_B);

p5e:[p5a,p5b,p5c,p5d];

result : ' charge balance ', above.
===========================

p6:subst(p5e,p4);

multiplication :V_Z+V_B
-------------------------------------

p6a:p6*(V_Z+V_B);

p6b:ratsimp(p6a);

set  'm:3*α(3)+2*α(2)+α(1)'
---------------------------------------

p6c:3*α(3)+2*α(2)+α(1) = m;

p6c1:rhs(p6c) = lhs(p6c);

p6d:subst([p6c],p6b);

find : formation factor  Φ: 'definition'
===============================

p7:Φ =(C_B*V_B)/(C_Z*V_Z);

p7a:lhs(p6d) - rhs(p6d)=0;

1e) find coef : V_Z,V_B
-----------------------------------

p7b :V_Z*ratcoef(p7a,V_Z) +V_B*ratcoef(p7a,V_B) ;

2e) solve V_B
---------------------

p7c: solve(p7b,V_B)[1]/V_Z;

4e) both sides : C_B/C_Z
---------------------------------------

p7d:p7c*(C_B/C_Z);

p7e:lratsubst([ rhs(p7)=lhs(p7)],(p7d));

insert  a: H-OH
----------------------

p7f:H-OH = a;

use : K_w =(OH)*(H)
------------------------------

p7f1:subst([p3a],p7f);

p7f1a:rhs(p7f1) = lhs(p7f1);

p7f2:lratsubst([ p7f],p7e);

num ,denom : ' see below'

p7f2a:ratsimp(num(rhs(p7f2))/(C_B*C_Z));
p7f2b:ratsimp(denom(rhs(p7f2))/(C_B*C_Z));

p7f2aa:expand(p7f2a);
p72bb:expand(p7f2b);

p8:lhs(p7f2) = p7f2aa/p72bb;

use : 'p6c1' and 'p7f1a'
-----------------------------------

p8a:subst([p6c1,p7f1a],p8);

use : definition , α(i)
------------------------------

p8b:subst([p2a,p2b,p2c,p2d],p8a);

how to find  :  d(d Φ/dph)/dph 'second order derivative ' 
================================================
find : inflection points  ' relation with EP '

d^2(Φ)/dph^2
============

1e)  ph = -log10(H)  => H = 10^(-ph)  , into  'p8b'
---------------------------------------------------------------------

p9:H = 10^(-ph);

p9a:lhs(p8b) = subst([p9],rhs(p8b));

2e) differentation
---------------------------

p9b:'diff(Φ,ph,2) = diff(rhs(p9a),ph,2)$

3e) back subst 'p9b'    ' 10^(-ph)  => H '
-------------------------------

p9c:lhs(p9b) =lratsubst([10^(-ph) = H],rhs(p9b))$

inflection points :d^2(Φ)/dph^2 = 0
=============================
rem :num ( )  = 0

p10:factor(num(rhs(p9c)));

values : 'see above' (H3PO4,NaOH).
==============================
H3PO4 :K_a1 : 10^(-2.23)
                K_a2 : 10^(-7.21)
                K_a3 :10^(-12.32)
                K_w : 10^(-14.00)

p10a:[C_B = 0.1,C_Z = 0.1,K_a1 =10^(-2.23),K_a2 =10^(-7.21),K_a3 =10^(-12.32),  K_w =10^(-14.00)];

p10b:subst(p10a,p10),numer,float;

-0.5301898110478399*H = 0  , '  H=0' , not physical solution
=================================================

We want search positive zeros of this polynomial
--------------------------------------------------------------------------

p11a:part(p10b,1,3);

roots
--------

p11a1:allroots(p11a);

p11a2:[H=2.202941921727431*10^-10,H=6.165921152838174*10^-8,H=1.800457439137481*10^-5,H=0.009875815640766335];

p11a3:[rhs(p11a2[1]),rhs(p11a2[2]),rhs(p11a2[3]),rhs(p11a2[4])];

1e) use : list ' all roots', pos physical meaning ' 

find corresponding :  H  to  ph = -log10(H)
---------------------------------------------------------------

p11aa:-log10(p11a3[1]),numer;
p11bb:-log10(p11a3[2]),numer;
p11cc:-log10(p11a3[3]),numer;
p11dd:-log10(p11a3[4]),numer;

rem
===

'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

real roots (positive)
--------------

p11b:realroots(p11a);

p11c:[p11b[4],p11b[5],p11b[6],p11b[7]],numer;

p11c1:[rhs(p11c[1]),rhs(p11c[2]),rhs(p11c[3]),rhs(p11c[4])];

find corresponding :  H  to  ph = -log10(H)
---------------------------------------------------------------

p12a:-log10(p11c1[1]),numer;
p12b:-log10(p11c1[2]),numer;
p12c:-log10(p11c1[3]),numer;
p12d:-log10(p11c1[4]),numer;

'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

ph : for 'inflection points'
=====================

p13:[p11aa,p11bb,p11cc,p11dd];

find : corresponding Φ  list for ph values
==================================
use : 'p9a','p10a'

p13a: subst(p10a,p9a)$

p13b: rhs(p13a)$
p13c: ΦΦ(ph) :=''p13b;

l1:[]$

for i:1 thru length(p13) step 1 do(
l1:endcons([ΦΦ(p13[i]),p13[i]],l1))$

first         coordinate : Φ value
=========================
second  coordinate : ph
=========================

l1;

mak graph :  ph = ph(Φ);
====================

ph : list
======

l2:makelist(i,i,0,13,0.1)$

l3:[];

for i:1 thru length(l2) step 1 do(
l3:endcons([ΦΦ(l2[i]),l2[i]],l3));

l3;

wxplot2d([[discrete, l3],[discrete,l3],[discrete,[l1[1]]],[discrete,[l1[2]]],[discrete,[l1[3]]],[discrete,[l1[4]]]]
    , [x ,-1,6], [y,0,14],
 [plot_format, xmaxima],[title,"titration weak acid,H3A '0.1M' with NaOH,'0.1M'"],
 [xlabel,"Φ"] , [ylabel,"ph"],
    [style,[points,0.7,red],[lines,0.6,blue],[points,4,blue],[points,4,blue],[points,4,blue],[points,4,blue]] , 
 [nticks,15],grid2d)$

see : blue mark on the graph are the inflections points  ? EP
==================================================
A = PO4, EP = Equivalence Point

relations
========

H3A      =  H2A(-) + H(+)     'K_a1'  = >  pK_a1 = 2.23
H2A(-)  =  HA(2-) + H(+)     'K_a2'  = >  pK_a2 =7.21
HA(2-)  =  A(2-)    + H(+)     'K_a3'  = >  pK_a3 =12.32

blue mark ( = points)
=================

p15:l1;

rem : Φ = 1  first  EP'1'
==================
1.000004597083253  =>  4.744617140422219  
rem : Φ = 2  first  EP'2'
==================
1.999969725206975  =>  9.656996952420625
rem : Φ = 3  first  EP'3'  ( approximation ,from list : l3)
===================
3.077607519109377  => 12.19999999999997 (assymptote)

p14:p10a;

pK_ai , i=1,2,3
============

p14a:-log10(rhs(p14[3])),float,numer;
p14b:-log10(rhs(p14[4])),float,numer;
p14c:-log10(rhs(p14[5])),float,numer;

p14d:[pK_a1:p14a,pK_a2:p14b,pK_a3:p14c];

approximation: EPi ,i=1,2,3 => (Φ = 1,2,3)
==================================
see : 'p15'

p16a:EP1=1/2*(pK_a1+pK_a2);  /*buffer approximation */
p16b:EP2=1/2*(pK_a2+pK_a3); /*buffer  approximation */
p16c:EP3=pK_a3;

'----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------'