Dear mathematical acolytes,
I am working as a materials scientist and my current topic is related to some probabilistic considerations of the microstructures of metals. I have the probability of something occuring $ P_\mathrm{occur} $ which is dependent on the diameter, $ c $, of something in the microstructure. This occurance probability is then modified with the probability density distribution of $c$ as $ P_\mathrm{occur}\left( c \right) p_\mathrm{size}\left( c \right) $ where
$$
P_\mathrm{occur}\left( c \right) = 1 – \exp{\left( -B\left( \frac{c}{c_\mathrm{N}} \right)^3 \right)},
\\
p_\mathrm{size}\left( c \right) = \frac{m-1}{c_\mathrm{min}}\left( \frac{c}{c_\mathrm{min}} \right)^{-m}.
$$
The probability density $ p_\mathrm{size}\left( c \right) $ is a power-law distribution with a lower bound at $ c_\mathrm{min} $ (a positive number) and upper bound at $ \infty $. The exponent is subject to $ m > 1 $. The parameter $ B $ in the occurance probability assumes a value between $ 0 $ and $ 1 $, and $c_\mathrm{N} $ is a positive number of the same order of $ c_\mathrm{min} $. The modified occurance probability is then to be integrated over the size $ c $ as:
$$
\int_{c_\mathrm{min}}^\infty{P_\mathrm{occur}\left( c \right) p_\mathrm{size}\left( c \right) \mathrm{d}c} = \\
\int_{c_\mathrm{min}}^\infty{ \left( 1 – \exp{\left( -B\left( \frac{c}{c_\mathrm{N}} \right)^3 \right)} \right) \frac{m-1}{c_\mathrm{min}}\left( \frac{c}{c_\mathrm{min}} \right)^{-m} {d}c}
$$
Here the first question appears, how is this to be integrated? Can it be shown using common mathematical functions?
As I am of little mathematical prowess, I tried using the symbolical library SymPy in Python to carry out the integration, and also to compare with a numerical integration using the trapezoidal rule. What happens is that the symbolical and the numerical integrations yields the same result, which is great, but that the symbolical integration returns a function named lowergamma(a,x), SymPy Docs, with a negative argument a evaluated at x. This is supposed to be the incomplete lower gamma Wikipedia, however, the incomplete lower gamma should not take negative arguments.
What puzzles me is that the symbolical integration with inserted numerical values in SymPy gives the same result as the trapezoidal numerical integration, even though the symbolical result uses the lower incomplete gamma with negative values, which should not be possible.
Here a cluster of questions appears, how can SymPy use negative values in its lower incomplete gamma? Neither SciPy nor Matlab can do this. What is different? How can I account for this in another library/language?
Thankful for comments and answers!
Update after answer from K.defaoite
The negative values in the argument of the incomplete gamma function also appears in the accepted answer below. This can be computed using 8.5.1 in the DLMF where a confluent hypergeometric function is used. Using this along with the solution presented below gives the same answer as both the symbolical and the numerical integration. The code is updated.
The code I used for my symbolical/numerical comparison outputs this:
Symbolic integral before evaluation = Symbolic integral before evaluation = 0.279982455532194*gamma(2/3)*lowergamma(-2/3, 0.5)/gamma(5/3) + 1.0 - 0.419973683298291*gamma(-2/3)
Value of symbolical integration = 0.744656471805332
Value of numerical integration = 0.7446564718064296
Integral by K.defaoite = 0.7446564718053317
The code itself, including the parameters used for evaluation:
import sympy as sym
import numpy as np
import scipy.special
# ---- Symbolical integration ----
c = sym.symbols('c', positive=True)
cmin = sym.symbols('cmin', positive=True)
cmax = sym.symbols('cmax', positive=True)
cN = sym.symbols('cN', positive=True)
B = sym.symbols('B', positive=True)
m = sym.symbols('m', positive=True)
func = (m-1)/(cmin**(1-m))*c**-m*(1 - sym.exp(-B*(c/cN)**3))
integrated = sym.integrate(func,(c, cmin, sym.oo), conds='separate')
# ---- Evaluation and comparison to numerical integration ----
Bnum = 0.5 # Somewhere between 0 and 1
cminnum = 0.1 # Small number, lower bound of power law dist.
cNnum = 0.1 # Parameter on the order of cmin
mnum = 3 # Exponent of power law dist, must be larger than 1, realistic up to 20
func = integrated[0].subs(B, Bnum).subs(cN, cNnum).subs(cmin, cminnum).subs(m, mnum)
cnum = np.logspace(np.log10(cminnum), 82, 100000000)
func_num = (mnum-1)/(cminnum**(1-mnum))*cnum**-mnum*(1. - np.exp(-Bnum*(cnum/cNnum)**3.))
trapz_func_num = np.trapz(func_num, cnum)
# ---- Integral by K.defaoite ----
a = (mnum - 1.)/3.
z = Bnum*(cminnum/cNnum)**3.
upper_gamma = scipy.special.gamma(-a) - z**(-a)/(-a)*scipy.special.hyp1f1(-a, -a+1., -z, out=None)
int_Kdefaoite = 1. - a*z**a*(upper_gamma)
# ---- Compare ----
print('Symbolic integral before evaluation =', func)
print('Value of symbolical integration = ', sym.N(func))
print('Value of numerical integration = ', trapz_func_num)
print('Integral by K.defaoite = ', int_Kdefaoite)
Best Answer
Make the change of variable $t=c/c_{\min}$ and let $z=B(c_{\min}/c_{\mathrm N})^3$ Which changes your integral to $$(m-1)\int_1^\infty \big(1-\exp(-z t^3)\big)t^{-m}\mathrm dt \\ =1+(1-m)\int_1^\infty \exp(-zt^3)t^{-m}\mathrm dt$$ Changing variables once more to $s=zt^3$ transforms this to $$1+(1-m)\int_z^\infty\exp(-s)\left(\frac{s}{z}\right)^{-m/3} \frac{\mathrm ds}{3s^{2/3}z^{1/3}} \\ =1+\frac{1-m}{3}z^{(m-1)/3}\int_z^\infty \exp(-s)s^{-(m+2)/3}\\ =1+\frac{1-m}{3}z^{(m-1)/3}\Gamma\left(\frac{1-m}{3},z\right) \\ =\boxed{1-az^{a}\Gamma\left(-a,z\right)}$$ Where $\Gamma(\cdot,\cdot)$ is the upper incomplete Gamma function and $a=(m-1)/3$ is positive. Unfortunately computing incomplete Gamma functions with negative arguments is quite hard - see the following paper and in particular reference 6.
EDIT: Some work has been done on asymptotic expressions for $\gamma(a,z)$ for large negative $a$ and $z$, but unfortunately, these expressions are not very accurate when $a$ is large and negative and $z$ is small and positive, as is the case in your problem.