Construct an array spacing proportional to a function or other array

Question

I have a function (f : black line) which varies sharply in a specific, small region (derivative f' : blue line, and second derivative f'' : red line). I would like to integrate this function numerically, and if I distribution points evenly (in log-space) I end up with fairly large errors in the sharply varying region (near 2E15 in the plot).

How can I construct an array spacing such that it is very well sampled in the area where the second derivative is large (i.e. a sampling frequency proportional to the second derivative)?

I happen to be using python, but I'm interested in a general algorithm.

---

Edit:
1) It would be nice to be able to still control the number of sampling points (at least roughly).
2) I've considered constructing a probability distribution function shaped like the second derivative and drawing randomly from that --- but I think this will offer poor convergence, and in general, it seems like a more deterministic approach should be feasible.

Answer 1

Assuming f'' is a NumPy array, you could do the following

# Scale these deltas as you see fit
deltas = 1/f''
domain = deltas.cumsum()

To account only for order of magnitude swings, this could be adjusted as follows...

deltas = 1/(-np.log10(1/f''))

Answer 2

I'm just spitballing here ... (as I don't have time to try this out for real)...

Your data looks (roughly) linear on a log-log plot (at least, each segment seems to be... So, I might consider doing a sort-of integration in log-space.

log_x = log(x)
log_y = log(y)

Now, for each of your points, you can get the slope (and intercept) in log-log space:

rise = np.diff(log_y)
run = np.diff(log_x)
slopes = rise / run

And, similarly, the the intercept can be calculated:

# y = mx + b
# :. b = y - mx
intercepts = y_log[:-1] - slopes * x_log[:-1]

Alright, now we have a bunch of (straight) lines in log-log space. But, a straight line in log-log space, corresponds to y = log(intercept)*x^slope in real space. We can integrate that easily enough: y = a/(k+1) x ^ (k+1), so...

def _eval_log_log_integrate(a, k, x):
    return np.log(a)/(k+1) * x ** (k+1)

def log_log_integrate(a, k, x1, x2):
    return _eval_log_log_integrate(a, k, x2) - _eval_log_log_integrate(a, k, x1)

partial_integrals = []
for a, k, x_lower, x_upper in zip(intercepts, slopes, x[:-1], x[1:]):
    partial_integrals.append(log_log_integrate(a, k, x_lower, x_upper))

total_integral = sum(partial_integrals)

^{You'll want to check my math -- It's been a while since I've done this sort of thing :-)}