Hi,
On Mar 11, 2011, at 2:17 AM, g-e-i-r-e wrote:
> Hi,
>
> I've been testing sympy for a phd project I'm working on and I've run
> into some challenges that I need some advice on. Consider the first
> part in the example code below, where the goal is to mix a function
> and some calculated parameter, utilizing a thermodynamic package which
> calculates some derivatives. First I define some constants/
> expressions which are necessary for my calculations and then some
> implicit functions representing what I need from the thermodynamic
> package (to make sure that I get the derivatives correctly). Finally,
> I combine them into the expression rNH3. Numerical values of this
> expressions and its derivatives with respect to T,V,n are needed for
> the calculations in my main routine (reactor calculations).
> In "CODE PART 2" I've given an example of how I calculate the
> quantity rNH3 and substitute in the values from the constants in "CODE
> PART 1" in addition to the ones from the thermodynamic package. In
> "CODE PART 3" I calculate the derivative of rNH3 with respect to T and
> substitute in the values of the constants and the thermodynamic
> package.
>
>
> Questions:
> 1) Can the code below be simplified in a manner that I still achieve
> my goal?
> 2) In "CODE PART 3" I need to make sure that I substitute in the
> expressions of thermodynamic package derivatives before I substitute
> in the non-derivatives of the thermodynamic package. Doing this in the
> opposite sequence will result in wrong results. This is probably due
> to ignorance, but does the subs function work as a pure textual search
> and replace of the expressions?
Well, technically no, because sometimes subs is a little mathematically smart,
like
In [74]: print (x**2).subs(x**3, y)
y**(2/3)
But other than that, it pretty much is. You can provide a dict or a list of
tuples to subs and it will perform all the substitutions you list. A list of
tuples is generally better because they will be performed in the order you give
them, where as a dict might be iterated through in any order. For example:
In [75]: (x*y).subs({x:z, y:7})
Out[75]: 7⋅z
In [76]: (x*y).subs([(x, z), (y, 7)])
Out[76]: 7⋅z
> 3) If there an easier path for substituting in values into the
> expressions? It would be nice to somehow tell a_NH3, a_N2, a_H2 their
> numerical values of the expressions and derivatives. Example:
> a_NH3 = 1.0 #
> a_NH3.diff(T) = 3.4 #
> a_NH3.diff(V) = 1.0 #
>
> and when the evalf function is called these values should be picked
> up. Comments?
You can give a substitution dictionary to evalf() that will be applied when
evaluating the expression. This is (I think) more efficient than substituting
before calling evalf().
In [80]: (x*y).evalf(subs={x:pi, y:7})
Out[80]: 21.9911485751286
Note that this only works with numerical values, or expressions that can be
evaluated to numerical values.
Aaron Meurer
>
> Thanks in advance!
>
>
> ----- CODE PART 1 -------
> rk0 = 8.849e14 # Rate constant for
> Arrhenius equation [kmol/h]
> kJprkcal = 4.1868 # Conversion factor -
> kJ / kcal
> rE = 40765.0 * kJprkcal # kJ / kmol
> rR = 8.314 * 1e3 # kJ / kmol kK
> rc = 0.5 # Constant -
> varies from 0.5 - 0.75
> k,k0,E,R,T,V,c = var('k k0 E R T V c')
> n = var('n')
> k = k0 * exp(-E / R / T) # Arrhenius equation
>
> a_NH3 = Function('a_NH3')(T, V, n) # Activity of
> ammonia. Retrieve from thermodynamic package
> a_N2 = Function('a_N2') (T, V, n) # Activity of
> Nitrogen. Retrieve from thermodynamic package
> a_H2 = Function('a_H2') (T, V, n) # Activity of
> Hydrogen. Retrieve from thermodynamic package
> Ka = Function('Ka')(T) # Equilibrium
> constant. Retrieve from thermodynamic package
>
> rNH3 = 2*k*(Ka**2.0*a_N2*(a_H2**3.0/a_NH3**2.0)**c -
> (a_NH3**2.0/a_H2**3.0)**(1.0-c))
>
>
> ---- CODE PART 2 -----
>
> def rate(t):
> substconsts = [(k0, rk0),(E, rE), (R, rR), (c, rc)]
> a = getactivity()
> Ka_ = getKaFunctionPtr()
> rNH3 = rNH3.subs([(a_NH3, a[0]),(a_N2, a[1]), (a_H2,
> a[2])])
> rNH3 = rNH3.subs([(Ka, Ka_(t))])
> rNH3 = rNH3.subs([(T, t)])
> rNH3 = rNH3.subs(substconsts)
> return 1e3*rNH3.evalf()
>
> -- CODE PART 3 ----
>
> def rate_x(t):
> a = getactivity()
> a_x = getactivity_x()
> Ka_ = getKaFunctionPtr()
> Ka_x = getKa_xFunctionPtr()
> substconsts = [(k0, rk0),(E, rE), (R, rR), (c, rc)]
>
> rNH3_x = rNH3.diff(T)
> rNH3_x = rNH3_x.subs([(a_NH3.diff(T), a_x[0]),
> (a_N2.diff(T), a_x[1]), (a_H2.diff(T), a_x[2])])
> rNH3_x = rNH3_x.subs([(Ka.diff(T), Ka_x(t))])
> rNH3_x = rNH3_x.subs([(a_NH3, a[0]),(a_N2, a[1]), (a_H2,
> a[2])])
> rNH3_x = rNH3_x.subs([(Ka, Ka_(t))])
> rNH3_x = rNH3_x.subs([(T, t)])
> rNH3_x = rNH3_x.subs(substconsts)
>
> return 1e3 * rNH3_x.evalf()
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to [email protected].
> To unsubscribe from this group, send email to
> [email protected].
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
