I tried it, and it was faster than the previous genetic algorithm. It took two minutes to calculate 500 points.
It is also possible that it takes time to calculate convex hull between algorithms.
The new method is based on a physical model, which makes n points on the sphere repel each other.
The repulsive force decreases with the increase of the distance between points, for example, it can be a force that decreases with the square of the distance.
This problem can be compared to the motion of n dots with equal charges on a sphere.
When the movement finally reaches balance and the balls stop moving, the distribution of the balls will be uniform due to the repulsive force.
At this time, the problem becomes a numerical problem of multivariable differential equations, and the final equilibrium state is what we want.
With the simplest Euler difference formula, six state variables of X, Y, Z, vx, vy and vz are set for each point to describe the motion state.
At each step, the motion state of a point is calculated, and the initial state is randomly assigned, with a speed of 0.
Function? []=main()
n = 12; ? Percentage of points
a=rand(N, 1)* 2 * pi; % is evenly distributed according to the random surface, and Mr enters the initial state.
b=asin(rand(N, 1)* 2- 1);
r0=[cos(a)。 *cos(b),sin(a)。 *cos(b),sin(b)];
V0 = zero (size (r0)););
g = 1e-2; % repulsion constant, this value is better.
For what? Ii= 1:200% simulates 200 steps, which generally converges, but in fact you can quit under it.
[rn,vn]=countnext(r0,v0,G); % update status
r0 = rnv0 = vn
end
' plot3(rn(:, 1),rn(:,2),rn(:,3),'.'); Keep? Open; % drawing results
[xx, yy, ZZ]= sphere (50);
h2=surf(xx,yy,ZZ); ? % Draw a unit ball for reference.
set(h2,' edgecolor ',' none ',' facecolor ',' r ',' facealpha ',0.7);
Axis? Equality;
Axis ([- 1? 1? - 1? 1? - 1? 1]);
Keep? Close;
end
Function? 【rn? vn]=countnext(r,v,G)? % Function to update status
%r stores the x, y and z data of each point, and v stores the speed data of each point.
num=size(r, 1);
Dd = zero (3, num, num); ? Percentage of vector difference between points
For what? M= 1: number-1
For what? N=m+ 1: quantity
dd(:,m,n)=(r(m,)-r(n,)';
dd(:,n,m)=-dd(:,m,n);
end
end
L=sqrt(sum(dd。 ^2, 1)); Percentage of distance between points
L(L & lt; 1e-2)= 1e-2; ? % distance too small limit
F = sum (dd. /repmat(L.^3,[3? 1? 1]),3)'; % calculated resultant force
Fr=r.*repmat(dot(F,r,2),[ 1? 3]); ? % Calculate the radial component of resultant force.
Fv = F-Fr; ? % tangential component
rn = r+v; % update coordinates
rn=rn。 /repmat(sqrt(sum(rn。 ^2,2)),[ 1? 3]);
VN = v+ G * Fv; % and new speed
end
Another advantage of this method is that you can see the results of each step.
The following is a process of finding the uniform distribution of 12 points by an algorithm.
The final icosahedron is very close to the regular icosahedron.