Optimizing spacing of mesh containing a given set of points

I tried to summarize the this as best as possible in the title. I am writing an initial value problem solver in the most general way possible. I start with an arbitrary number of initial values at arbitrary locations (inside a boundary.) The first part of my program creates a mesh/grid (I am not sure which is the correct nuance), with N points total, that contains all the initial values. My goal is to optimize the mesh such that the spacing is as uniform as possibl开发者_开发百科e. My solver seems to work half decently (it needs some more obscure debugging that is not relevant here.)

I am starting with one dimension. I intend to generalize the algorithm to an arbitrary number of dimensions once I get it working consistently. I am writing my code in fortran, but feel free to reply with pseudocode or the language of your choice.

Allow me to elaborate with an example:

Say I am working on a closed interval [1,10]

xmin=1  
xmax=10 

Say I have 3 initial points: xmin, 5 and xmax

num_ivc=3  
known(num_ivc)=[xmin,5,xmax] //my arrays start at 1. Assume "known" starts sorted

I store my mesh/grid points in an array called coord. Say I want 10 points total in my mesh/grid.

N=10  
coord(10)  

Remember, all this is arbitrary--except the variable names of course. The algorithm should set coord to {1,2,3,4,5,6,7,8,9,10}

Now for a less trivial example:

num_ivc=3  
known(num_ivc)=[xmin,5.5,xmax  
or just  
num_ivc=1  
known(num_ivc)=[5.5]   

Now, would you have 5 evenly spaced points on the interval [1, 5.5] and 5 evenly spaced points on the interval (5.5, 10]? But there is more space between 1 and 5.5 than between 5.5 and 10. So would you have 6 points on [1, 5.5] followed by 4 on (5.5 to 10]. The key is to minimize the difference in spacing.

I have been working on this for 2 days straight and I can assure you it is a lot trickier than it sounds. I have written code that

only works if N is large
only works if N is small
only works if it the known points are close together
only works if it the known points are far apart
only works if at least one of the known points is near a boundary
only works if none of the known points are near a boundary

So as you can see, I have coded the gamut of almost-solutions. I cannot figure out a way to get it to perform equally well in all possible scenarios (that is, create the optimum spacing.)

Your initial points define a number of intervals whose lengths are (say) a1,a2,...,ak. You're going to divide them into (say) m1,m2,...,mk subintervals. Clearly you want equal spacing within each interval, so you'll have m1 subintervals of length a1/m1, then m2 of length a2/m2, etc.

Now, of course you haven't defined "optimum spacing" and there probably isn't any One True Way to define it, but let's take it to mean "minimizing the sum of the squares of the subinterval lengths". (Easy exercise: when you have only one interval, this does in fact give you an equal subdivision.) Then your aim is to minimize a1^2/m1 + ... + ak^2/mk, with the constraint that m1+...+mk is fixed and the m's are all positive integers.

Unfortunately, discrete optimization is hard. What's the answer if you let the m's vary continuously? Unsurprisingly it's that the m's should be proportional to the a's. Here's a simple approach that should do tolerably well (I'm guessing you don't strictly need an optimal solution).

1. Compute an initial guess by simply rounding the "continuous" solution to the nearest integers. That is: if you're aiming for a total of M+1 points, hence M subintervals, take mj = round(M*aj/s) where s = a1+...+ak.
2. Unfortunately this won't give you the correct value of m1+...+mk unless you happen to be lucky. So now fix that up by either adding 1 to some of them (if you've got too few points) or subtracting 1 from some of them (if you've got too many). Do this, say, in order of the error incurred by rounding; that is, in order of M*aj/s - round(M*aj/s). (Increasing or decreasing order depending on which direction your m1+...+mk was wrong in.)

There's an assumption I'm making here that may or may not be correct, namely that what you want is something like "all intervals of equal length, as nearly as possible" rather than "adjacent intervals of equal length, as nearly as possible". (Test case: suppose you're trying to subdivide the interval [0,1] with 9 internal points and you have to have one of them at 0.001. Do you take 0, 0.001, and 8 more points equally spaced between 0.001 and 1, or do you actually want them more closely spaced near 0.001?) If it's the latter you want, then things get more interesting and more difficult.

And a word of warning: Do not expect anything you do in one dimension to generalize neatly to give something that works well in more dimensions. In dimensions higher than 1, you'll probably want not just a set of points but something like a triangulation, and (1) that has a lot more structure than your one-dimensional setup and (2) the conditions you'll want on the triangles are more complicated -- e.g., you probably won't care only about the size; you'll also want them not to be too long and thin.