Build an array interval proportional to a function or another array

Question

I have a function ( f : black line) that changes dramatically in a specific small area (derivative f' : blue line and second derivative f'' : red line). I would like to integrate this function numerically, and if I distribute the points evenly (in the log space), I get quite large errors in a sharply changing area (about 2E15 on the graph).

How can I construct an array interval so that it is very well selected in the region where the second derivative is large (i.e., the sampling frequency proportional to the second derivative)?

I am using python, but I'm interested in the general algorithm.

---

Edit:
1) It would be nice to be able to control the number of sample points (at least approximately).
2) I considered the question of building a probability distribution function formed as a second derivative, and drawing randomly from this --- but I think this will offer poor convergence, and in general it seems a more deterministic approach should be feasible.

Answer 1

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

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

To take into account only fluctuations of the order of magnitude, this could be changed as follows:

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

Answer 2

I'm just spitballing here ... (since I don’t have time to actually try this) ...

Your data looks (roughly) linear on the log chart (at least each segment seems ... So, I could consider integrating into the log space.

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

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

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

And, similarly, the interception can be calculated:

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

Ok, now we have a bunch of (direct) lines in the log-log space. But the straight line in the log-log space corresponds to y = log(intercept)*x^slope in real space. We can easily integrate this: 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 want to check my math - some time has passed since I did such things :-)}

Answer 3

1) Cool approach

At the moment, I have adopted the “adaptive refinement” approach inspired by hydrodynamic methods . I have a function that I want to try, f , and I select some initial array of sample points x_i . I create a "fetch" g function that determines where to insert new fetch points.

In this case, I chose g as the slope of log(f) ---, since I want to allow quick changes in the log space . Then I divide the range g into levels L=3 refinement. If g(x_i) exceeds the refinement level, this range is subdivided into N=2 pieces, these units are added to the samples and checked to the next level. This gives something like this:

The solid gray line is the function that I want to sample, and the black crosses are my starting points of the sample.
The dashed gray line is a derivative of the log of my function.
colored dashed lines are my "refinement levels"
colored crosses are my specified sample points.

All this is displayed in the log space.

2) Simple approach

After I finished (1), I realized that maybe I just chose the maximum distance in y and choose x -spacings to achieve this. Similarly, just divide the function evenly by y and find the corresponding points x . The results of this are shown below:

Answer 4

A simple approach would be to split the x-axis array into three parts and use different intervals for each of them. This will allow you to maintain the total number of points, as well as the necessary interval in different regions of the site. For instance:

 x = np.linspace(10**13, 10**15, 100) x = np.append(x, np.linspace(10**15, 10**16, 100)) x = np.append(x, np.linspace(10**16, 10**18, 100)) 

You can choose the best interval based on your data, but you will get this idea.