pH of 40% solution of sodium bisulfite in water (init. pH=7)

Simplistic_Eq_NaHSO3.py

from sympy.matrices import Matrix
from sympy import solve, nsolve, symbols, pprint
import numpy as np
from numpy import log10
import timeit
from Tkinter import Text, Tk, END, INSERT

# i			C0, gmol/L
# 1	C0_{H2SO3} = 	0
# 2	C0_{HSO3(-)} =	5.10E+00
# 3	C0_{SO3(2-)} =	0
# 4	C0_{H3O(+)} = 	1.00E-07
# 5	C0_{OH(-)} =	1.01E-07
# 6	C0_{H2O} =	5.54E+01
# 7	C0_{Na(+)} = 	5.10E+00

# r1	H2SO3	+H2O<==>H3O(+)	+HSO3(-)	Kc1*Ce6 = 10^-pKa1 = Ce4*Ce2/Ce1=10^-1.8569852
# r2	HSO3(-)	+H2O<==>H3O(+)	+SO3(2-)	Kc2*Ce6 = 10^-pKa2 = Ce4*Ce3/Ce1=10^-7.171984936
# r3	H2O	+H2O<==>OH(-)	+H3O(+)		Kc3*Ce6^2 = 10^-pKw = Ce4*Ce5=10^-13.99567863
#								(Harris)
# 	Ci = Ci0 + Sum(nu_ij csi_j)											

# 	ECS					VARS
# 1	0 =-Ce1 -xe1				Ce1
# 2	0 =-Ce2 + C02 + xe1 - xe2		Ce2
# 3	0 =-Ce3 + xe2				Ce3
# 4	0 =-Ce4 + C04 + xe1 + xe2 + xe3		Ce4
# 5	0 =-Ce5 + C05 + xe3			Ce5
# 6 	0 =-Ce6 + C06 - xe1 - xe2 -2xe3		Ce6
# 7 	0 =-Ce7 + C07				Ce7
# 6	0 = Ka1*Ce1 - Ce4*Ce2			xe1
# 7	0 = Ka2*Ce2 - Ce4*Ce3			xe2
# 8	0 = Kw - Ce4*Ce5			xe3

def main():
	eps = np.finfo(float).eps
	C1,C2,C3,C4,C5,C6,C7,x1,x2,x3 = symbols('C1,C2,C3,C4,C5,C6,C7,x1,x2,x3')
	X0 = (0.01,5.08,0.01,1e-7,10**-13.99567863/1e-7,997/(2*1.00794+15.9994),5.10,0,0,0)
	X = (C1,C2,C3,C4,C5,C6,C7,x1,x2,x3)	
	f_0 = (
		-C1-x1,
		-C2+5.100858754+x1-x2,
		-C3+x2,
		-C4+1e-7+x1+x2+x3,
		-C5+10**-13.99567863/1e-7+x3,
		-C6 -x1 -x2 -2*x3 + 997/(2*1.00794+15.9994),
		-C7+5.100858754,
		-10**-1.8569852/(997/(2*1.00794+15.9994))+C4*C2/(C1*C6),
		-10**-7.171984936/(997/(2*1.00794+15.9994))+C4*C3/(C2*C6),
		-10**-13.99567863/(997/(2*1.00794+15.9994))**2+C4*C5/C6**2)
	J = Matrix((
		(-1, 0, 0, 0, 0, 0, 0,-1, 0, 0),
		(0, -1, 0, 0, 0, 0, 0,+1,-1, 0),
		(0, 0, -1, 0, 0, 0, 0, 0,+1, 0),
		(0, 0, 0, -1, 0, 0, 0,+1,+1,+1),
		(0, 0, 0, 0, -1, 0, 0, 0, 0,+1),
		(0, 0, 0, 0, 0, -1, 0,-1,-1,-2),
		(0, 0, 0, 0, 0,  0,-1, 0, 0, 0),
		(-C2*C4/(C1**2*C6),C4/(C1*C6),0,C2/(C1*C6),0,-C4*C2/(C1*C6**2),0,0,0,0),
		(0,-C3*C4/(C2**2*C6),C4/(C2*C6),C3/(C2*C6),0,-C4*C3/(C2*C6**2),0,0,0,0),
		(0,0,0,C5/C6**2,C4/C6**2,-2*C4*C5/C6**3,0,0,0,0))
		)
	pprint(J)
	print '\n'
	nu_ij = np.matrix([[-1, 0, 0], [1, -1, 0], [0, 1, 0], [1, 1, 1], [0, 0, 1], [-1, -1, -2], [0, 0, 0]])
	print '\n'
	J=Matrix(np.zeros([10,10]).astype(int))
	J[0:7, 0:7] = -1*np.eye(7).astype(int)
	J[0:7, 7:10] = nu_ij
	J[7:10, 0:7] = ( \
		np.diag(np.power(np.matrix(X[0:7]).T, -1).A1,0) * \
		nu_ij * \
		np.diag(np.prod(np.power(np.matrix(X[0:7]).T, nu_ij), 0).A1) \
		).T
	text2 = pprint(J)	
	solution = nsolve(f_0, X, X0)
	a 		= np.empty(5,dtype='S20')
	a[0] 	= 'pH =' + '{:.6f}'.format(-log10(float(solution[3])))
	a[1] 	= 'pOH =' + '{:.6f}'.format(-log10(float(solution[4])))
	a[2] 	= 'pH + pOH =' + '{:.6f}'.format(-log10(float(solution[3]))+-log10(float(solution[4])))	
	a[3] 	= 'CH2O0 =' + '{:.6f}'.format(float(X0[5]))
	a[4] 	= 'CH2O =' + '{:.6f}'.format(float(solution[5]))
	root = Tk()
	T = Text(root, height=40, width=30)
	T.pack()
	for part in a:
		print part
		T.insert(END, str(part) + '\n')
	T.insert(END, 'Solution ' + str(X) + '= \n')
	for part in solution:
		T.insert(END, str(part)+', \n')
	T.delete('end -1 lines',END)
	T.delete('end -3 chars',END)
	T.delete('end -3 chars',END)
	T.insert(END, '\n'*3)
	X0 = [0.0111594604012,\
5.07850930681,\
0.0111899867912,\
3.05436900904e-05,\
1.83001028601e-08,\
55.3418793956,\
5.100858754,\
-0.0111594604012,\
0.0111899867912,\
-8.26998962602e-08]
	T.insert(END, 'X0 = ' + str(X0))
	T.insert(END, '\n'*2+'f(X0) =')
	replacements = [(X[i],X0[i]) for i in range(10)]
	sum_f_2 = 0
	subs_part = 0
	for part in f_0:
		subs_part = part.subs(replacements)
		sum_f_2 += subs_part**2
		T.insert(END, '\n'+str(subs_part))
	T.insert(END, '\n'*1+'||f(X0)|| =' + str(sum_f_2))

main()
raw_input("Press Enter to continue...")