GreenGroup/RMG-Py — Issue #62: Numerical issues with Wilhoit-toMultiNASA conversion, apparently associated with use of indefinite analytic integrals
Author: gmagoon
Created: 2011-11-12 03:33:41+00:00
State: open

I've come across a tricky case that seems to cause numerical problems in the Wilhoit-to-MultiNASA conversion:

>>> w2 = thermo.Wilhoit(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> ans = thermo.Wilhoit_integral_T0(w2,6.0) - thermo.Wilhoit_integral_T0(w2,.298)
>>> print ans
-27.3531739832
>>> print (thermo.Wilhoit_integral_TM1(w2,1.5) -thermo.Wilhoit_integral_TM1(w2,.298))
69.7868424873
>>> print (thermo.Wilhoit_integral2_TM1(w2,1.5) -thermo.Wilhoit_integral2_TM1(w2,.298))
3263.51989762
>>> print thermo.TintOpt_objFun_W(0.65, w2, .298, 1.5, 3)
-8.5055398813e-07
>>> print thermo.TintOpt_objFun_W(0.60, w2, .298, 1.5, 3)
-1.57702925208e-06
>>> print thermo.TintOpt_objFun_W(0.6037, w2, .298, 1.5, 3)
2.38989377976e-06
>>> print thermo.TintOpt_objFun_W(0.615, w2, .298, 1.5, 3)
8.49388015922e-07
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, True, 3)
ERROR:root:Greg thought he fixed the numerical problem, but apparently it is still an issue; please e-mail him with the following results:
ERROR:root:0.6
ERROR:root:Wilhoit(cp0=4, cpInf=75, a0=-54.497423456539401, a1=303.4090028535806, a2=-727.01542657068092, a3=659.08598716992981, H0=1234.5, S0=789, B=3)
ERROR:root:0.298
ERROR:root:1.5
ERROR:root:True
ERROR:root:-1.57702925208e-06
0
>>> print thermo.Wilhoit_integral2_TM1(w2,1.5)
-574582.356995
>>> print thermo.Wilhoit_integral2_TM1(w2,.298)
-577845.876893

The above shows:
*Negative values of the objective function with fairly significant magnitude; these are inaccurate as the objective function is an integral of squared error
*Jittery, non-smooth behavior of the objective function
*Significant magnitude of the indefinite integrals, which are used to compute the definite integrals (by difference)

My investigation of this issue suggests is that it ultimately arises from numerical issues associated with the use of intermediate indefinite analytic integrals with large magnitude in the calculation. It may be possible to implement definite analytic integrals that avoid differencing of large magnitude intermediate values.

Using an older version of RMG-Py that has numerical integral capabilities, I have confirmed that use of definite numerical integrals smooths the objective function and keeps it positive.

Using the older version with the numerical option (definite numerical integrals in both objective function and NASA coefficient calculation):

>>> print thermo.Cp_TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
8.57655311464e-005
>>> print thermo.Cp_TintOpt_objFun(0.6037, w2, .298, 1.5, 1, 3)
8.38349646131e-005
>>> print thermo.Cp_TintOpt_objFun(0.615, w2, .298, 1.5, 1, 3)
7.93030143242e-005
>>> print thermo.Cp_TintOpt_objFun(0.630, w2, .298, 1.5, 1, 3)
7.61815561034e-005
>>> print thermo.Cp_TintOpt_objFun(0.640, w2, .298, 1.5, 1, 3)
7.57683296015e-005
>>> print thermo.Cp_TintOpt_objFun(0.650, w2, .298, 1.5, 1, 3)
7.65980407849e-005
>>> print thermo.Cp_TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.00101810388993
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.000430559256529

Using a modified version of the aforementioned which uses definite numerical integrals for just the objective function (indefinite analytic integrals used in matrix for NASA coefficient calculation):

>>> from rmg import thermo
>>> w2 = thermo.ThermoWilhoitData(4.0,75.0,-54.497423456539401,303.4090028535806,-727.01542657068092,659.08598716992981,1234.5,789.0,B=3.0)
>>> print (w2.integral_T0(6.0) - w2.integral_T0(.298))
-27.3531739832
>>> print thermo.TintOpt_objFun(0.60, w2, .298, 1.5, 1, 3)
7.56744975661e-006
>>> print thermo.TintOpt_objFun(1.0, w2, .298, 1.5, 1, 3)
0.00043583716888
>>> print thermo.TintOpt_objFun(0.61, w2, .298, 1.5, 1, 3)
7.57669477025e-006
>>> print thermo.TintOpt_objFun(0.62, w2, .298, 1.5, 1, 3)
7.7480199252e-006
>>> print thermo.TintOpt_objFun(0.63, w2, .298, 1.5, 1, 3)
8.12499638414e-006
>>> print thermo.TintOpt_objFun(0.64, w2, .298, 1.5, 1, 3)
8.66185655468e-006
>>> print thermo.TintOpt_objFun(0.59, w2, .298, 1.5, 1, 3)
7.68816880736e-006
>>> print thermo.TintOpt_objFun(0.58, w2, .298, 1.5, 1, 3)
8.1157850218e-006
>>> print thermo.TintOpt_objFun(0.603, w2, .298, 1.5, 1, 3)
7.51797415433e-006
>>> print thermo.TintOpt_objFun(0.607, w2, .298, 1.5, 1, 3)
7.57045654609e-006

A problem with similar symptoms was encountered in issue #20, but this was a cython issue; this issue does not appear to be related to cython, and instead seems to be a more fundamental numerical difficulty.