With[{n = {3, 4, 5, 6, 7, 20}},
Partition[show[points[#]] & /@ n, 3]] // GraphicsGrid
With the code below, with piecewise spring forces (exclusively repulsive) and Cartesian vectors. A different model for force and switch to spherical vectors would probably be wise.
show[points_] :=
Graphics3D[{
Opacity[.5], Sphere[],
Opacity[1], PointSize[Medium], Point[points]},
ImageSize -> Tiny,
Boxed -> False]
points[n_] :=
With[{
steps = 50,
start = RandomPoint[Sphere[], n]},
Nest[move[#, .1] &, start, steps]]
move[points_, dt_] :=
With[{n = Length[points]},
Module[{copy = points},
Do[
copy[[i]] += dt Sum[force[copy[[j]], copy[[i]]],
{j, DeleteCases[Range[n], i]}];
copy[[i]] = copy[[i]]/Norm[copy[[i]]],
{i, n}];
copy]]
force[p1_, p2_] :=
With[{d = VectorAngle[p1, p2], range = Pi},
If[d > range, 0, Normalize[p2 - p1]*(range - d)]]
Update: Simulated Annealing
Start with n random points on the unit sphere. Then do the main loop: pick a point randomly and make a random move. If the move lowers the energy (electro-static potential), accept the move into a new configuration. Accept the move also, this is crucial, if the energy difference is less then Exp[-difference/T], where T is a temperature-like control parameter. Do K1 * Length[points] of these inner iterations (about 10, so each point is picked about ten times) and K2 outer passes of the inner iterating. Before each outer pass T is set, i.e. gradually lowered from some starting value (few times higher than a typical temperature change) to zero.
Results agree with literature on this, so-called Thomson problem.
start = points[18];
result = Reap[simulation[start, 10, 10]]; // AbsoluteTiming
(* {2.7106, Null} *)
Energy is recorded in the main loop:
ListPlot[result[[2, 1]],
Joined -> True,
ImageSize -> Medium]
Table[Graphics3D[{
PointSize[Medium], Point[xyz /@ p],
Opacity[.75], Sphere[{0, 0, 0}, .99]},
ImageSize -> Small,
Boxed -> False,
SphericalRegion -> True],
{p, {start, result[[1]]}}] // GraphicsRow
Code:
points[n_Integer] :=
With[{x := RandomReal[]},
Table[{2 Pi \[Xi], ArcCos[2 \[Xi] - 1]}, n]]
xyz[{u_, v_}] :=
{Sin[v] Cos[u], Sin[v] Sin[u], Cos[v]}
energy[points_] :=
With[{n = Length[points]},
Sum[1/EuclideanDistance @@ (xyz /@ points[[{i, j}]]),
{i, n - 1}, {j, i + 1, n}]]
(* O(n) instead of energy O(n^2) *)
dE[points_, i_, r_] :=
With[{n = Length[points]},
Module[{indices, before, after},
indices = DeleteCases[Range[n], i];
before = Sum[
1/Norm[xyz@points[[j]] - xyz@points[[i]]], {j, indices}];
after = Sum[
1/Norm[xyz@points[[j]] - xyz@r], {j, indices}];
after - before]]
simulation[start_, K1_, K2_] :=
With[{
n = Length[start],
du = .1 Pi,
dv = .1 Pi,
x := RandomReal[]},
Module[{T, temp = start},
Do[
T = 10 (1 - (k1 - 1)/K1)^5;
temp = Nest[
Module[{i, p, change},
(* monitor energy; must Reap *)
Sow[energy[#]];
i = RandomInteger[{1, n}];
p = #[[i]] + {du (2 x - 1), dv (2 x - 1)};
change = dE[#, i, p];
If[change < 0 || (x < Exp[-change/T]),
ReplacePart[#, i -> p], #]] &, temp, K2 n], {k1, K1}];
temp]]
