Skip to content

Instantly share code, notes, and snippets.

@santiago-salas-v

santiago-salas-v/test_davies_2.ipynb

Created July 9, 2016 14:40
Show Gist options
  • Save santiago-salas-v/39543056f74267dd32a253bbea6720f8 to your computer and use it in GitHub Desktop.
Save santiago-salas-v/39543056f74267dd32a253bbea6720f8 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
test_davies_2.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<pre>\n",
"```\n",
"| i | comp. | m0, mol/g |\n",
"|---|-------|-----------|\n",
"|0|H2O\t\t|\t55.5084350618\n",
"|1|H2SO3\t|\t0.0\n",
"|2|HSO3(-)\t|\t0.0\n",
"|3|SO3(2-)\t|\t0.0\n",
"|4|H3O(+)\t|\t1.00300902708e-10\n",
"|5|HO(-)\t|\t1.01303910853e-10\n",
"|6|Na(+)\t|\t5.11620737613\n",
"```\n",
"</pre>\n",
"$$x = [n_0 m_1 m_2 m_3 m_4 m_5 m_6 \\xi_1 \\xi_2 \\xi_3 \\gamma_0 \\gamma_1 \\gamma_3 \\gamma_4 \\gamma_5 \\gamma_6]$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from sympy import *\n",
"import numpy as np\n",
"from IPython.display import display\n",
"init_printing(use_latex='mathjax')\n",
"eps = np.finfo(float).eps\n",
"mm0 = 2*1.00794+15.9994\n",
"m0_ref = 1/1000.0\n",
"x, n0, m1, m2, m3, m4, m5, m6 = symbols('x n0 m1 m2 m3 m4 m5 m6')\n",
"xi1, xi2, xi3 = symbols('xi1 xi2 xi3')\n",
"g0, g1, g2, g3, g4, g5, g6 = symbols('g0 g1 g2 g3 g4 g5 g6')\n",
"i = symbols('i')\n",
"# X0 = ( 55.3418794136,\n",
"# 1.1193032265e-05,\n",
"# 0.00509379348046,\n",
"# 1.12236685315e-05,\n",
"# 3.06355951606e-08,\n",
"# 3.31668819804e-13,\n",
"# 0.00511621018126,\n",
"# -0.0111594470496,\n",
"# 0.0111899913906,\n",
"# -1.00669325488e-07,\n",
"# 0.306009594112,\n",
"# 0.92620164755,\n",
"# 0.735905750795,\n",
"# 0.92620164755,\n",
"# 0.92620164755,\n",
"# 1.00118133956,\n",
"# 0.92620164755,\n",
"# 0.00512746448589)\n",
"X0 = (5.53419097566066e+01,\n",
" 2.22712743154495e-26,\n",
" 5.11620737612838e-03,\n",
" 2.22712743154495e-26,\n",
" 1.00300902708124e-10,\n",
" 1.01303910852863e-10,\n",
" 5.11620737612838e-03,\n",
" 0.00000000000000e+00,\n",
" 0.00000000000000e+00,\n",
" 0.00000000000000e+00,\n",
" 8.31653143040483e-01,\n",
" 1.00117874448073e+00,\n",
" 9.26272473142106e-01,\n",
" 7.36130872136064e-01,\n",
" 9.26272473142106e-01,\n",
" 9.26272473142106e-01,\n",
" 9.26272473142106e-01,\n",
" 5.11620747693079e-03)\n",
"X = (n0, m1, m2, m3, m4, m5, m6,\n",
" xi1, xi2, xi3,\n",
" g0, g1, g2, g3, g4, g5, g6, i)\n",
"Q1 = 1 * g0**-1 * g1**-1 * g2**1 * g3**0 * g4**1 * g5**0 * g6**0 *\\\n",
" (1/mm0)**-1 * m1**-1 * m2**1 * m3**0 * m4**1 * m5**0 * m6**0 *\\\n",
" (1/m0_ref)**(-1-1+1+0+1+0+0)\n",
"Q2 = 1 * g0**-1 * g1**0 * g2**-1 * g3**1 * g4**1 * g5**0 * g6**0 *\\\n",
" (1/mm0)**-1 * m1**0 * m2**-1 * m3**1 * m4**1 * m5**0 * m6**0 *\\\n",
" (1/m0_ref)**(-1+0-1+1+1+0+0)\n",
"Q3 = 1 * g0**-2 * g1**0 * g2**0 * g3**0 * g4**1 * g5**1 * g6**0 *\\\n",
" (1/mm0)**-2 * m1**0 * m2**0 * m3**0 * m4**1 * m5**1 * m6**0 *\\\n",
" (1/m0_ref)**(-2+0+0+0+1+1+0)\n",
"f_0 = (\n",
" -n0 + 55.3418794136 - 1*xi1 - 1*xi2 -2*xi3,\n",
" -m1*n0*mm0 + 2.22712743154495e-23/1000.0*(55.3418794136)*mm0 - xi1,\n",
" -m2*n0*mm0 + 5.11620737612839e+00/1000.0*(55.3418794136)*mm0 + xi1 - xi2,\n",
" -m3*n0*mm0 + 2.22712743154495e-23/1000.0*(55.3418794136)*mm0 + xi2,\n",
" -m4*n0*mm0 + 1.00300902708124e-07/1000.0*(55.3418794136)*mm0 + xi1 + xi2 + xi3,\n",
" -m5*n0*mm0 + 1.01303910852863e-07/1000.0*(55.3418794136)*mm0 + xi3,\n",
" -m6*n0*mm0 + 5.11620737612839e+00/1000.0*(55.3418794136)*mm0,\n",
" -10**-1.8569852/(1/(mm0 * m0_ref))**-1 + Q1,\n",
" -10**-7.171984936/(1/(mm0 * m0_ref))**-1 + Q2,\n",
" -10**-13.99567863/(1/(mm0 * m0_ref))**-2 + Q2,\n",
" -g0 + exp(-1.0*mm0*(m1+ m2 + m3 + m4 + m5 + m6)),\n",
" -g1 + 10**(0.1*i/m0_ref),\n",
" -g2 + 10**(-0.510*(-1.0)**2*(sqrt(i/m0_ref)/(1+sqrt(i/m0_ref))-0.3*i/m0_ref)),\n",
" -g3 + 10**(-0.510*(-2.0)**2*(sqrt(i/m0_ref)/(1+sqrt(i/m0_ref))-0.3*i/m0_ref)),\n",
" -g4 + 10**(-0.510*(+1.0)**2*(sqrt(i/m0_ref)/(1+sqrt(i/m0_ref))-0.3*i/m0_ref)),\n",
" -g5 + 10**(-0.510*(-1.0)**2*(sqrt(i/m0_ref)/(1+sqrt(i/m0_ref))-0.3*i/m0_ref)),\n",
" -g6 + 10**(-0.510*(+1.0)**2*(sqrt(i/m0_ref)/(1+sqrt(i/m0_ref))-0.3*i/m0_ref)),\n",
" -i + 1/2.0*(m2*(-1)**2 + m3*(-2.0)**2 + m4*(+1.0)**2 + m5*(-1.0)**2 + m6*(+1.0)**2)\n",
" )\n",
"solution = nsolve(f_0, X, X0, verbose=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Matrix([[-1, 0, 0, 0, 0, 0, 0, -1, -1, -2, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m1, -18.01528*n0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m2, 0, -18.01528*n0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m3, 0, 0, -18.01528*n0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m4, 0, 0, 0, -18.01528*n0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m5, 0, 0, 0, 0, -18.01528*n0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [-18.01528*m6, 0, 0, 0, 0, 0, -18.01528*n0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], \n",
" [0, -18.01528*g2*g4*m2*m4/(g0*g1*m1**2), 18.01528*g2*g4*m4/(g0*g1*m1), 0, \n",
" 18.01528*g2*g4*m2/(g0*g1*m1), 0, 0, 0, 0, 0, -18.01528*g2*g4*m2*m4/(g0**2*g1*m1), \n",
" -18.01528*g2*g4*m2*m4/(g0*g1**2*m1), 18.01528*g4*m2*m4/(g0*g1*m1), 0, \n",
" 18.01528*g2*m2*m4/(g0*g1*m1), 0, 0, 0], \n",
" [0, 0, -18.01528*g3*g4*m3*m4/(g0*g2*m2**2), 18.01528*g3*g4*m4/(g0*g2*m2), \n",
" 18.01528*g3*g4*m3/(g0*g2*m2), 0, 0, 0, 0, 0, -18.01528*g3*g4*m3*m4/(g0**2*g2*m2), \n",
" 0, -18.01528*g3*g4*m3*m4/(g0*g2**2*m2), 18.01528*g4*m3*m4/(g0*g2*m2), \n",
" 18.01528*g3*m3*m4/(g0*g2*m2), 0, 0, 0], \n",
" [0, 0, -18.01528*g3*g4*m3*m4/(g0*g2*m2**2), 18.01528*g3*g4*m4/(g0*g2*m2), \n",
" 18.01528*g3*g4*m3/(g0*g2*m2), 0, 0, 0, 0, 0, -18.01528*g3*g4*m3*m4/(g0**2*g2*m2), \n",
" 0, -18.01528*g3*g4*m3*m4/(g0*g2**2*m2), 18.01528*g4*m3*m4/(g0*g2*m2), \n",
" 18.01528*g3*m3*m4/(g0*g2*m2), 0, 0, 0], \n",
" [0, -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" -18.01528*exp(-18.01528*m1 - 18.01528*m2 - 18.01528*m3 - 18.01528*m4 - 18.01528*m5 - 18.01528*m6), \n",
" 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 100.0*10**(100.0*i)*log(10)], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, \n",
" 10**(-16.1276160668587*sqrt(i)/(31.6227766016838*sqrt(i) + 1) + 153.0*i)* \\\n",
" (153.0 + 255.0/(31.6227766016838*sqrt(i) + 1)**2 - 8.06380803342937/ \\\n",
" (sqrt(i)*(31.6227766016838*sqrt(i) + 1)))*log(10)], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, \n",
" 10**(-64.5104642674349*sqrt(i)/(31.6227766016838*sqrt(i) + 1) + 612.0*i)* \\\n",
" (612.0 + 1020.0/(31.6227766016838*sqrt(i) + 1)**2 - 32.2552321337175/ \\\n",
" (sqrt(i)*(31.6227766016838*sqrt(i) + 1)))*log(10)], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, \n",
" 10**(-16.1276160668587*sqrt(i)/(31.6227766016838*sqrt(i) + 1) + 153.0*i)* \\\n",
" (153.0 + 255.0/(31.6227766016838*sqrt(i) + 1)**2 - 8.06380803342937/ \\\n",
" (sqrt(i)*(31.6227766016838*sqrt(i) + 1)))*log(10)], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, \n",
" 10**(-16.1276160668587*sqrt(i)/(31.6227766016838*sqrt(i) + 1) + 153.0*i) \\\n",
" *(153.0 + 255.0/(31.6227766016838*sqrt(i) + 1)**2 - 8.06380803342937/ \\\n",
" (sqrt(i)*(31.6227766016838*sqrt(i) + 1)))*log(10)], \n",
" [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, \n",
" 10**(-16.1276160668587*sqrt(i)/(31.6227766016838*sqrt(i) + 1) + 153.0*i)* \\\n",
" (153.0 + 255.0/(31.6227766016838*sqrt(i) + 1)**2 - 8.06380803342937/ \\\n",
" (sqrt(i)*(31.6227766016838*sqrt(i) + 1)))*log(10)], \n",
" [0, 0, 0.500000000000000, 2.00000000000000, 0.500000000000000, 0.500000000000000,\n",
" 0.500000000000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1]])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment