9
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ I found this link (my.originlab.com/forum/topic.asp?TOPIC_ID=1251) which seems to head the same way I would normally go: compute slope and intercept in log-log space. Then convert to lin-lin space, integrate, and evaluate. $\endgroup$ – MrMas Jun 21 '18 at 15:17
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for your response! But I think this is effectively just performing the integral in linear space. Perhaps I'm asking for something impossible however... $\endgroup$ – DilithiumMatrix Oct 2 '15 at 3:26
  • 2
    $\begingroup$ No, this does the integral in log space. When discretizing, $da$ is equally sized in log space, not linear space. $\endgroup$ – Damascus Steel Oct 2 '15 at 4:45
  • 1
    $\begingroup$ @DilithiumMatrix is right: The discretisation of the $x$ values are in log-space, but the interpolation of the $y$ values happens in linear space. Thus, if you were to use the trapezoidal rule, the function that is effectively integrated is piecewise linear in a plot with logarithmic x-axis and linear y-axis. $\endgroup$ – burnpanck Apr 30 '18 at 15:09
3
$\begingroup$

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$.)

$\endgroup$
  • $\begingroup$ Thanks for your (very clear) response, but as I've just said in response to @DamascusSteel --- I think this is just reverting the integral to linear-linear space, and losing the benefits of log-space. $\endgroup$ – DilithiumMatrix Oct 2 '15 at 3:27
  • 1
    $\begingroup$ @DilithiumMatrix: This is not the same answer as the one by DamascusSteel. Note that applying trapezoid-rule to DamascusSteel's answer would not give the $\frac{\exp(b)-1}{b}$ factor. $\endgroup$ – burnpanck Apr 30 '18 at 15:18
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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: enter image description here $$\log(y_1)=m_1 \log(x_1) + n_1$$ $$\log(y_2)=m_1 \log(x_1) + 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.

$\endgroup$
  • $\begingroup$ This is wonderful @elenapascal, this has been bothering me for over 3 years now, and I think this is (or is very close to) the solution. I don't quite following your last relation, I don't think the integral over y equal to the log(x2/x1) $\endgroup$ – DilithiumMatrix Apr 4 at 22:53
  • $\begingroup$ In particular, if I take the log of the integral on the left-hand-side, then I get a similar term to the right hand side, but with log([x_2/x_1]^{m_1+1} + 1), i.e. there is an additional +1 in the argument of the log $\endgroup$ – DilithiumMatrix Apr 4 at 23:03
  • $\begingroup$ It bothered me a lot today as well, only after I wrote it I realised @Stefan had posted the same answer. For m=-1 you just replace that in the definition of y: y(x)=e^n/x. That gives logs. I'm not sure I follow your second post $\endgroup$ – Elena Pascal Apr 4 at 23:09
  • $\begingroup$ I just realized the same thing, but I hadn't fully understood until I read your explanation $\endgroup$ – DilithiumMatrix Apr 4 at 23:11
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.