I am currently using build_function from Symbolics.jl to generate Julia-type functions that I later work with (numerically integrate, etc.). Recently, I was not able to numerically integrate one of such-obtained functions and I soon realized, why. Fixing all but one variables the function in question has the following form (I write it the same way Julia returns it to me):
f(\varepsilon) = b_1 + \dfrac{b_2}{\left[1 + \dfrac{k_1 \left( \varepsilon^{2} \right)^{0.375}}{ \exp\left(c_1/\varepsilon^{2}\right) - 1} \right]^2 \left[1 + \dfrac{k_2 \left( \varepsilon^{2} \right)^{0.375}}{\exp\left(c_2/\varepsilon^{2}\right) - 1} \right]^2} \\ + b_3 \varepsilon \dfrac{\left( 2 + \dfrac{k_3 \left( \varepsilon^{2} \right)^{0.375}}{\exp\left(c_3/\varepsilon^{2}\right) - 1}\right) \left( \dfrac{g_1 \varepsilon}{\left[\exp\left(c_3/\varepsilon^{2}\right) - 1\right] \left( \varepsilon^{2} \right)^{0.625}} + \dfrac{g_2 \varepsilon \exp\left(c_4/\varepsilon^{2}\right)}{\left[\exp\left(c_4/\varepsilon^{2}\right) - 1\right]^{2} \left( \varepsilon^{2} \right)^{1.625}} \right)}{\left\lbrace \left[1 + \dfrac{k_5 \left( \varepsilon^{2} \right)^{0.375}}{\exp\left(c_5/\varepsilon^{2}\right) - 1} \right]^2\right\rbrace^{2} \left[1 + \dfrac{k_6 \left( \varepsilon^{2} \right)^{0.375}}{\exp\left(c_6/\varepsilon^{2}\right) - 1} \right]^2}
with a_i, b_i, and c_i (\simeq 1) being (positive) constants, some of which might be equal to each other. In any case, it’s rather evident that \lim_{\varepsilon\to0} f(\varepsilon) = b_1 + b_2, which in this case \simeq 0.03. This is nicely illustrated by the plot below:
It also, however, shows that Julia struggles to evaluate f(\varepsilon) when \varepsilon gets too small: in the dashed region, f returns NaN. The problem seems to lie in the expression
\dfrac{g_2 \varepsilon \exp\left(c_4/\varepsilon^{2}\right)}{\left[\exp\left(c_4/\varepsilon^{2}\right) - 1\right]^{2} \left( \varepsilon^{2} \right)^{1.625}}
since in the \varepsilon\to0 limit the latter has the Inf/Inf form (even though it’s clear that the denominator’s Inf is the “stronger” one).
Is there any workaround? The problem would probably be solved if I could factor out the leading-order exponential in each term. Perhaps, Symbolics.jl has some built-in function that deals with such limits, which I was missed in the documentation…
P.S. Perhaps, minimal (non-)working example could help:
using Symbolics
@variables x;
f(x) = x*exp(1/x^2)/(x^3.25*(exp(1/x^2) - 1)^2);
g = build_function(f(x),x,expression=Val{false});
g(0.03)