Like the title says, I'm interested in numerically evaluating the function $$I(x)=\frac{\int_0^xf(x^\prime)\mathrm{d}x^\prime}{x}$$ where $f(x^\prime)$ is given by the user as sampled points and is greater than zero at $x=0$ so that this function evaluates to a finite value. The trouble is that near $x=0$ both the numerator and denominator become small and lead to trouble, especially when you are actually at $x=0$.

I could obviously Taylor expand in this region, but using a piecewise expression feels like a kluge and means I can't fully take advantage of vectorized math to speed things up due to the if-statement.

Is there some global way of re-writing this function that makes it safe for numerical evaluation near $x=0$?

  • $\begingroup$ Taylor expansion of $F=\int f$ should give good approximations. The first order approximation gives $f(0)$ $\endgroup$ – Maximilian Janisch Feb 3 at 16:47
  • $\begingroup$ Also, is Taylor expansion not „vectorizable“ (what is the programming context to your question?) $\endgroup$ – Maximilian Janisch Feb 3 at 16:49
  • $\begingroup$ Could you expand on "$f(x′)$ is given by the user as sampled points"? So it means you have only the values $f(k\Delta x)$ available, and would for example integrate the quadratic interpolation polynomial using the first three points for the first two intervals? $\endgroup$ – Lutz Lehmann Feb 3 at 16:50
  • $\begingroup$ @MaximilianJanisch, by not vectorizable I mean that if I Taylor expand then my code would look like def I(x): if(x<epsilon): do_something(x) else: do_something_else(x) The problem is that I want to evaluate $I(x)$ between $10^4$ and $10^6$ times and care about performance for my application so the if statement means I can't take advantage of vectorized operations acting over the array of $x$ in something like numpy. That's easily a factor of three or more in execution time. $\endgroup$ – ElectronsAndStuff Feb 3 at 17:08
  • $\begingroup$ @LutzLehmann, that's just a low-level detail that maybe I shouldn't have emphasized the way I did. It's exactly what you said where I have the values of the function at $x=k\Delta x$ for non-negative integers $k$. Just off the top of my head I would be using numpy's implementation of the trapezoidal rule. $\endgroup$ – ElectronsAndStuff Feb 3 at 17:11

We can use the MVT for integrals to get that $$\int_0^x f(x')\mathrm{d}x'=xf(c)$$ For some $c \in (0,x)$. This guarentees that $$\lim_{x \to 0} I(x)=f(0)$$ If $f$ is continuous.

| cite | improve this answer | |
  • $\begingroup$ The assumption of continuity of $f$ is needed for the first statement (the mean-value theorem for integrals). $\endgroup$ – John Bentin Feb 3 at 16:38
  • $\begingroup$ Thanks! That's an interesting point. I guess the difficulty then lies in finding $c(x)$ since in general the point will change for each different interval you use. $\endgroup$ – ElectronsAndStuff Feb 3 at 17:00
  • $\begingroup$ @JohnBentin I wanted to express that continuity is necessary for the whole process with "If $f$ is continuous." at the end. But I'm going to put it to the front to make it clear. $\endgroup$ – Botond Feb 3 at 17:46

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.