Integral in log-log space

Question

I'm working with functions which, in general, are much smoother and better behaved in log-log space --- so that's where I perform interpolation/extrapolation, etc, and that works very well. Is there a way to integrate these numerical functions in log-log space?

i.e. I'm hoping to use some sort of simple trapezoidal rule to perform a cumulative integral (e.g. in python, use scipy.integrate.cumtrapz), to find some $F(r)$ s.t.

$$F(r) = \int_0^r y(x) \, dx$$

But I'm hoping to use the values $log(y)$ and $log(x)$, instead of $y$ and $x$ (when possible).

Answer 1

I think that there is a bit of confusion with change of variables in some of the previous answers as well as some errors. The integral of a log function is not the log of the integral. I think in general is difficult to write out the the integral of a function knowing the integral of its log. If anyone knows how to do that I would be interested.

In the meanwhile, @Stefan's solution above is the way to get around integrating a function in log-log space. The starting point is that the function you're dealing is linear in log-log space for small enough segments.

One can then write the equation of the line at the segment endpoints: $$\log(y_1)=m_1 \log(x_1) + n_1$$ $$\log(y_2)=m_1 \log(x_2) + n_1$$

where $m_1$ is the slope of the line and $n_1$ is its $y$-intercept.

By subtracting the two, one can find:

$$m_1=\frac{\log(y_1) -\log(y_2)}{\log(x_1)-\log(x_2)}$$

and from substitution: $$n_1=\log(y_1)-m_1 \log(x_1)$$

If in the log-log space the equation of a segment is close to a line then in normal (linear) space the equation of the segment is close to an exponential:

$$y(x)\simeq x^{m} e^{n}$$

If we have a analytical formulation for this segment it is easy to integrate:

$$\int_{x_1}^{x_2}y(x)dx = \frac{e^{n_1}}{m_1+1}(x_2^{m_1+1}-x_1^{m_1+1}), \,\, \text{for } \, m\neq -1$$

and $$\int_{x_1}^{x_2}y(x)dx = e^{n_1} \log \frac{x_2}{x_1},\, \,\text{for } \, m = -1 $$

This feels a bit like cheating, but this is sampling in log-log space such that we can approximate the function in the linear space to an exponential with parameters derived from the log-log space.

Answer 2

You can just change variables. Setting $a=log(x)$, $b(a)=log(y(x))$. The integral becomes

$F(r)=\int^{log(r)}_{-\infty} exp(a+b) da$

You have to be a little careful because you are integrating from $-\infty$. What you have to do exactly will depend on what $y(x)$ looks like.

Answer 3

I do not use python, but if I understand correctly then by $$ F(r)=\int_0^ry(x)dx $$ you are thinking something like $$ \mathbf{F}=\mathrm{integrate}(\mathbf{y},\mathbf{x}) $$ where $\mathbf{F}=[F_1,...,F_n]$ is a vector sampling the integral over a grid $\mathbf{x}$.

However you do not have samples of $x$ and $y$, but rather you have samples of $\hat{x}=\log(x)$ and $\hat{y}=\log(y)$.

Of course the simplest approach would be $$ \mathbf{F}=\mathrm{integrate}(\exp(\mathbf{\hat{y}}),\exp(\mathbf{\hat{x}})) , $$ but this would be error-prone, because $y(x)$ is not smooth, even though $\hat{y}(\hat{x})$ is.

Now the trapezoidal rule essentially assumes your input $y(x)$ is piecewise linear. So the simple generalization would be for you to assume that $\hat{y}(\hat{x})$ is piecewise linear.

In this case, defining $\Delta F_k=F_{k+1}-F_k$, you have $$ \Delta F_k=\int_{x_k}^{x_{k+1}}y(x)dx =\int_{\hat{x}_k}^{\hat{x}_{k+1}}e^{\hat{y}(\hat{x})}e^\hat{x}d\hat{x} =\int_{\hat{x}_k}^{\hat{x}_{k+1}}\tilde{y}(\hat{x})d\hat{x} $$

Then, defining $t=(\hat{x}-\hat{x}_k)/\Delta \hat{x}_k$, you have $$ \hat{y}_{k+t} \approx \hat{y}_k + t \Delta \hat{y}_k $$ and $\tilde{y}(t) \approx a e^{bt}$, with $a=e^{\hat{y}_k+\hat{x}_k}$ and $b=\Delta \hat{y}_k+\Delta \hat{x}_k$.

So the integral becomes $$ \Delta F_k \approx a \Delta \hat{x} \int_0^1e^{bt}dt = a \Delta \hat{x} \frac{e^b-1}{b} $$

In Matlab this would look something like

dlogx=diff(logx); dlogy=diff(logy); k=1:length(logx)-1;  
b=dlogx+dlogy; a=exp(logx+logy);  
dF=a(k).*dlogx.*(exp(b)-1)./b;  
F=cumsum([0,dF]);

Hope this helps!

(Edit: My answer is essentially the same as the much more concise answer that Damascus Steel gave as I was typing. The only difference is I attempted to give a particular solution for the case where the "particular $y(x)$" is a piecewise-linear $\hat{y}(\hat{x})$ discretized over a discrete $\mathbf{\hat{x}}$ mesh, with $F(\hat{x}_1)=0$.)

Answer 4

If the function looks smooth in a log-log plot, you can interpolate using a power law on each interval (power laws are, of course, linear in log-log). Thus, between points $(x_i, y_i)$ and $(x_{i+1}, y_{i+1})$ under the assumption that $y = C_i x^{n_i}$ within interval $i$, you obtain $n_i = \ln(y_i/y_{i+1})/\ln(x_i/x_{i+1})$ and $C_i = \ln(y_i) - n_i\ln(x_i)$. The contribution to the integral from interval $i$ is then $$\Delta F_i = \int_{x_i}^{x_{i+1}} C_i x^{n_i}dx = \left\{ \begin{array}{l@{\quad}l} \frac{C_i}{n_i + 1} (x_{i+1}^{n_i+1} - x_i^{n_i+1}), & n_i \neq -1 \\ C_i(\ln x_{i+1} - \ln x_i), & n_i=-1, \end{array} \right. $$ where you obviously need some tolerance for identifying the special case $n_i=-1$ in your implementation.

Answer 5

The solution I use is basically an implementation of the trapezium rule and makes use of the scipy.misc.logsumexp function to maintain precision. If you have some function lny that returns the logarithm of y then you can do, e.g.:

from scipy.misc import logsumexp
import numpy as np

xmin = 1e-15
xmax = 1e-5

# get values of x spaced logarithmically
xvs = np.logspace(np.log10(xmin), np.log10(xmax), 10000)

# evaluate your function at xvs
lys = lny(xvs)

# perform trapezium rule integration
deltas = np.log(np.diff(xvs))
logI = -np.log(2.) + logsumexp([logsumexp(lys[:-1]+deltas), logsumexp(lys[1:]+deltas)])

The value logI is the log of the integral that you want.

Obviously this will not work if you need to set xmin = 0. But, if you have some non-zero positive lower limit to the integral you can just play about with the number of points in xvs to find a number where the integral converges.