6
$\begingroup$

I am trying to solve an equation including a definite integral. I have a real valued function:

$f(a,b,x) = \frac{(1-x)^a}{x}\exp(-b/x)$

Given $a_1, b_1$ and $b_2$, I wish to find $a_2$ from the condition:

$\int_0^1 dx\ f(a_1,b_1,x) = \int_0^1 dx\ f(a_2,b_2,x)$.

If we are worried about undefined behavior at $x=0$, we can set $f(0) = 0$ always.

My best attempt so far is to rewrite the condition as:

$\int_0^1 dx\ \ln(f(a_1,b_1,x)) = \int_0^1 dx\ \ln(f(a_2,b_2,x))$, which gives me: $a_2 = a_1 + \int_0^1 dx\ \frac{b_1 - b_2}{x}$. That is divergent, so I will never get a number out of it.

I can of course solve it numerically, which I in fact already do. It sits in an optimization routine. But it would be much faster if an analytic answer is possible.

$\endgroup$
2
  • $\begingroup$ For the integral to converge at $x =1$, $a>-1$. The integral is a monotonically decreasing function for $a>0$ and monotonically increasing function for $-1<a<0$. So the binary search isn't fast enough? $\endgroup$
    – anecdote
    Sep 20, 2016 at 19:49
  • 2
    $\begingroup$ you cannot rewrite the condition by taking the ln of the integrand - this will give a very different result than the orginal condition - and hence the divergent behavior. $\endgroup$
    – Andreas
    Sep 21, 2016 at 20:26

2 Answers 2

4
$\begingroup$

Some insights...

We can look at your integral as a Mellin convolution.

The Mellin transform is defined as

$$ \mathcal{M}\left\{f\right\}(s) = \mathcal{M}_f(s) = \int_0^\infty f(x)\,x^{s-1}dx$$

and the Mellin convolution of two functions is defined as

$$ (f *_{\tiny \mathcal{M}} g)(s) = \int_0^{\infty} f(x)\,g\!\left(\frac{s}{x}\right)\frac{dx}{x}. $$

In your case, we have $f(x) = (1-x)^a \theta(1-x)$, $g(x) = e^{-x}$, and we want to compute $(f *_{\tiny \mathcal{M}} g)(b)$.


A nice property of Mellin convolution is

$$ (f *_{\tiny \mathcal{M}} g)(s) = \mathcal{M}^{-1} \left\{\mathcal{M}_f(s) \mathcal{M}_g(s) \right\}(s), $$

and using the table here, we have

$$ \mathcal{M}_f(s)\mathcal{M}_g(s) = \frac{\Gamma(a+1)\Gamma^2(s)}{\Gamma(s+a+1)}, \quad \text{for} \;\; \text{Re}(a) > -1, \; \text{Re}(s) > 0. $$

Now by definitions of the inverse Mellin transform and the Meijer G function, we get

$$ (f *_{\tiny \mathcal{M}} g)(b) = \Gamma (a+1) G_{1,2}^{2,0}\left(b\left| \begin{array}{c} a+1 \\ 0,0 \\ \end{array} \right.\right).$$


We can apply the second identity here to simplify the Meijer G and express your integral as

$$ \boxed{\int_0^1 dx \; f(a, b, x) = e^{-b} \Gamma (a+1) U(a+1,1,b),} $$

where $U$ is the confluent hypergeometric function

$$ U(\alpha, \beta, z) = \frac{1}{\Gamma(\alpha)} \int_0^{\infty } e^{-zt}t^{\alpha-1} (t+1)^{\beta-\alpha-1} dt. $$

Other than identities here for specific $a$ and $b$, I'm not aware of a way to whittle down $U(a+1,1,b)$ to anything nicer.


Edit

Upon looking at the definition of $U(\alpha, \beta, z)$, I realize one can make the substitution $x \mapsto 1/(t+1)$ to get

\begin{align} \int_0^1 (1-x)^a \exp\left(-\frac{b}{x}\right) \frac{dx}{x} &= \int_0^\infty e^{-b(t+1)} t^a (t+1)^{-a-1} dt\\ &= e^{-b^{\phantom{1}^\phantom{1}}} \!\!\! \Gamma (a+1) U(a+1,1,b), \end{align}

which is a much less roundabout solution.

$\endgroup$
2
  • $\begingroup$ Yes, I have actually not gone through the trouble to go over the Mellin transform myself, but got to somewhat similar result by trying to look up the integral and play around a bit to try and ease the pain of solving it numerically. Unfortunately it does not become much easier to solve in the end :-). But thank you very much for the effort, I appreciate it. $\endgroup$
    – nabla
    Sep 23, 2016 at 7:15
  • $\begingroup$ Yep, exactly. It seems you need a solution that's (easily) symbolically invertible with respect to $a$. I'd wager it doesn't exist... $\endgroup$
    – Greg Hurst
    Sep 23, 2016 at 14:22
3
+50
$\begingroup$

I don't believe there is going to be an analytic solution for solving the integral for a. (I'm happy if someone proves me wrong).

So we need to do something practical. Here is some consideration.

You want to solve for $a_2$:

$$\int_0^1 dx\ f(a_2,b_2,x) - \int_0^1 dx\ f(a_1,b_1,x) = 0$$

Let's rename $- \int_0^1 dx\ f(a_1,b_1,x) = G_{-1}$.

With some initial good starting value, replace $a_2$ with $a_2 + d$ where the difference $d$ will be small. So you need

$$ G(d) = G_{-1} + \int_0^1 dx\ f(a_2+d,b_2,x) = 0$$

Now $f(a_2+d,b_2,x) = f(a_2,b_2,x) (1-x)^d $ and you can expand $(1-x)^d$ into a power series in $d$. This works, since you expand, in the integrand,

$$ (1-x)^d = \sum_{n=0}^\infty (\ln(1-x) )^n d^n / n! $$

EDIT: there was a wrong extra factor $(1-x)^d$ in the summands which was removed here and later.

and the terms do not diverge in the interval $x \in (0, 1)$ and the corresponding integrals exist as well. So you have

$$ G(d) = G_{-1} + \sum_{n=0}^\infty d^n \quad 1/ n! \int_0^1 dx\ f(a_2,b_2,x) (\ln(1-x) )^n = 0$$

Exemplarily, I chose some $a_2,b_2$ and I let Wolfram Alpha compute the first $n=10$ many terms $1/ n! \int_0^1 \cdots$ and one sees their values decrease rapidly with $n$.

So this leaves you with choosing some $N$ (depending on the desired accuracy), computing the integrals and then solving, for $d$, a polynomial equation $$ G(d) = G_{-1} + \sum_{n=0}^N d^n G_n = 0$$

Replacing $a_2$ with $a_2 +d$, this can then be iterated.

Hope this makes life easier ...

$\endgroup$
1
  • $\begingroup$ Nice! I like that! I will implement it right now and see if it saves me a couple of precious iterations... $\endgroup$
    – nabla
    Sep 24, 2016 at 15:11

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .