# Numerically Evaluating Functions Like $I(x)=(\int_0^xf(x^\prime)\mathrm{d}x^\prime)/x$ near $x=0$

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$$?

• Taylor expansion of $F=\int f$ should give good approximations. The first order approximation gives $f(0)$ – Maximilian Janisch Feb 3 at 16:47
• Also, is Taylor expansion not „vectorizable“ (what is the programming context to your question?) – Maximilian Janisch Feb 3 at 16:49
• 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? – Lutz Lehmann Feb 3 at 16:50
• @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. – ElectronsAndStuff Feb 3 at 17:08
• @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. – 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.
• The assumption of continuity of $f$ is needed for the first statement (the mean-value theorem for integrals). – John Bentin Feb 3 at 16:38
• 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. – ElectronsAndStuff Feb 3 at 17:00
• @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. – Botond Feb 3 at 17:46