12
$\begingroup$

I'm trying to program a random walk on a sphere with Mathematica. I found this code while I was searching, but I'm not getting results from it. I waited an entire day, but my PC didn't finish evaluating rw. I would like some help please (:

rotateWithAxis[p_, a_, theta_] := 
      #/Norm[#] & @ ((1 – Cos[theta]) (a.p) a + p Cos[theta] + Cross[a, p] Sin[theta]);

Using the function, the random walk on a unit sphere is written as follows:

rw = With[{stepLength = 0.03, num = 10000},
         Module[{rotateWithAxis, p, a, q},
         rotateWithAxis[p_, a_, theta_] := 
         #/Norm[#] &@((1 – Cos[theta]) (a.p) a + p Cos[theta] + Cross[a, p] Sin[theta]);
         a = {1., 0, 0}; q = {Cos[stepLength], Sin[stepLength], 0};
         Table[p = a; a = q; q = rotateWithAxis[p, a, RandomReal[{0, 2 Pi}]], {num}]]];
$\endgroup$
  • 3
    $\begingroup$ I don't know where did you get this from but in (1-Cos[theta]) the minus is in fact \[Dash]. Just rewrite this expression and it should work. $\endgroup$ – Kuba Mar 22 '14 at 11:45
  • $\begingroup$ p.s. minor edit: #/Norm[#] & is Normalize. $\endgroup$ – Kuba Mar 22 '14 at 11:47
  • $\begingroup$ my minus is [Dash], when i write that mathematica replace it with - didnt work anyway :c but thanks (: $\endgroup$ – Mariana da Costa Mar 22 '14 at 12:47
  • $\begingroup$ I'm not quite sure how have you managed to copy the code here without revealing this? While posting an answer I had to manually change \[Dash] to -. $\endgroup$ – Kuba Mar 22 '14 at 13:12
  • $\begingroup$ I do not agree that it is simple mistake. :) $\endgroup$ – Kuba Mar 23 '14 at 8:06
8
$\begingroup$

As noted by previous answerers, the desired distributional properties of the spherical random walk was not properly clarified. Nevertheless, let me offer two variations of interest.

The first variation is the spherical analog of the bounded random walk (this recent thread shows a few ways on how to implement this). This would seem to have been the variation that was being attempted. Before I can show my solution, let me pull out a few auxiliary routines:

(* https://mathematica.stackexchange.com/a/10994 *)
arc[center_?VectorQ, {start_?VectorQ, end_?VectorQ}] := Module[{ang, co, r},
    ang = VectorAngle[start - center, end - center];
    co = Cos[ang/2]; r = EuclideanDistance[center, start];
    BSplineCurve[{start, center + r/co Normalize[(start + end)/2 - center], end}, 
                 SplineDegree -> 2, SplineKnots -> {0, 0, 0, 1, 1, 1},
                 SplineWeights -> {1, co, 1}]]

(* slightly faster than the equivalent RotationMatrix[{vv1, vv2}];
   from http://dx.doi.org/10.1080/10867651.1999.10487509 *)
vectorRotate[vv1_?VectorQ, vv2_?VectorQ] := 
 Module[{v1 = Normalize[vv1], v2 = Normalize[vv2], c, d, d1, d2, t1, t2},
        d = v1.v2;
        If[TrueQ[Chop[1 + d] == 0],
           c = UnitVector[3, First[Ordering[Abs[v1], 1]]];
           t1 = c - v1; t2 = c - v2; d1 = t1.t1; d2 = t2.t2;
           IdentityMatrix[3] - 2 (Outer[Times, t2, t2]/d2 - 
           2 t2.t1 Outer[Times, t2, t1]/(d2 d1) + Outer[Times, t1, t1]/d1),

           c = Cross[v1, v2];
           d IdentityMatrix[3] + Outer[Times, c, c]/(1 + d) - LeviCivitaTensor[3].c]]

Here is a function that takes a bounded random step in the sphere.

boundedRandomStep[v_?VectorQ, φ_?NumericQ] := vectorRotate[{0, 0, 1}, v].({0, 0, Cos[φ]} + 
       Sin[φ] Append[Normalize[RandomVariate[NormalDistribution[], 2]], 0])

φ here is the length of the arc connecting the unit vector v and the generated random variate.

From this, here is how one might generate a bounded random walk:

With[{start = {0, 0, 1}, steps = 200, φ = π/15}, 
     BlockRandom[SeedRandom[42, Method -> "Legacy"]; 
                 Graphics3D[{Sphere[], {Red, Sphere[start, Scaled[1/150]]},
                             {Directive[Blue, Arrowheads[Small]], 
                              Arrow[Tube[arc[{0, 0, 0}, #], Scaled[1/1000]]] & /@ 
                              Partition[NestList[boundedRandomStep[#, φ] &, start, steps],
                                        2, 1]}}, Boxed -> False, PlotRange -> 1.2]]]

bounded random walk on a sphere


Another variation rests on using a distribution biased towards a particular "mean direction"; one such distribution is the the von Mises-Fisher distribution. First, here is a routine for generating von Mises-Fisher variates (previously shown in this answer):

vonMisesFisherRandom[μ_?VectorQ, κ_?NumericQ] := Module[{ξ = RandomReal[], w},
        w = 1 + (Log[ξ] + Log[1 + (1 - ξ) Exp[-2 κ]/ξ])/κ;
        RotationTransform[{{0, 0, 1}, Normalize[μ]}][
        Append[Sqrt[1 - w^2] Normalize[RandomVariate[NormalDistribution[], 2]], w]]]

Here is the random walk based on von Mises-Fisher:

With[{start = {0, 0, 1}, steps = 200, κ = 8}, 
     BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; 
     Graphics3D[{Sphere[], {Red, Sphere[start, Scaled[1/150]]},
                 {Directive[Blue, Arrowheads[Small]], 
                  Arrow[Tube[arc[{0, 0, 0}, #], Scaled[1/1000]]] & /@ 
                  Partition[NestList[vonMisesFisherRandom[#, κ] &, start, steps],
                            2, 1]}}, Boxed -> False, PlotRange -> 1.2]]]

von Mises-Fisher random walk

$\endgroup$
21
$\begingroup$

Here is an approach.

rws[n_, p0_?(Norm@# == 1 &), ang_] := 
 NestList[RotationMatrix[ang/(2 Pi),
          (Function[{u, v}, {Cos[u] Cos[v], Cos[u] Sin[v], Sin[u]}] @@
           RandomReal[{0, 2 Pi}, 2])].# &, p0, n]

Visualizing:

Graphics3D[{Sphere[], Line[rws[10000, {1, 0, 0}, 1]]}, Boxed -> False]

spherical random walk from different views

trace of spherical random walk

$\endgroup$
  • 2
    $\begingroup$ I am not convinced that this is correct. You seem to be trying to rotate around a random axis by ang/(2Pi) (why the 2Pi division at all?). But this random axis doesn't have a uniform distribution on the surface of the sphere (it has a higher density close to the poles). $\endgroup$ – Szabolcs Mar 22 '14 at 15:34
  • 1
    $\begingroup$ Also, rotating around an axis randomly chosen from the surface doesn't seem to be what's required. The rotation should be around a random axis chosen uniformly from the great circle defined by the plan perpendicular to the vector corresponding to the current point. $\endgroup$ – Szabolcs Mar 22 '14 at 15:35
  • $\begingroup$ @Szabolcs thank you...I agree it depends sensitively on what one means by random walk...as a naive first approach I merely imagined myself as an ant on the surface of perfectly symmetric surface, no preferred direction/isotopic and walking in regular angular step sizes. The rotation of the axis just reflects random choice of direction., the sphere rotating randomly under me is equivalent to me randomly walking. If a particular distribution is aim or other definition I look forward to seeing $\endgroup$ – ubpdqn Mar 23 '14 at 1:31
  • $\begingroup$ @Szabolcs tha ang/2$\pi$ was unnecessary just picked step size of radians...could have left open...could also randomize starting point...thank you again and look forward to other and better approaches $\endgroup$ – ubpdqn Mar 23 '14 at 1:39
  • $\begingroup$ @Szabolcs ...and yes I appreciate that as a little I chose geodesic paths not the 'realistic' walks but it was a nice simplifying assumption...albeit unrealistic and with anomalous consequences $\endgroup$ – ubpdqn Mar 23 '14 at 1:47
12
$\begingroup$

If something is evaluating forever, it is because of poor implementation, an error or because it is just too much work to do :).

What I'm doing while debugging, is: if it looks ok -> run the minimal example, if it fails -> run it step by step:

screenshot

As you can see Abs[0. - 0.0987745-] is an expression that has no way to appear if - was a real Minus.

Fortunately, there aren't many explicitly written subtractions; you can copy each and try:

(1 – Cos[theta]) // FullForm
Times[\[Dash], Cos[theta]]
$\endgroup$
  • $\begingroup$ yes, i'm gonna do that thanks :) $\endgroup$ – Mariana da Costa Mar 22 '14 at 20:59
4
$\begingroup$

A neat compact solution for self-avoiding random walk on a sphere. No explicit coordinate transformation, - all based on built in algorithms.

pts=RandomPoint[s=Sphere[{0,0,0},1],1000];
walk=pts[[Last[FindShortestTour[pts]]]];

That's it. To visualize I will not use arcs of sphere but simple line segments, as only sequence and positions of points are important usually and lines look pretty close to arcs for a large number of points. Let me know, if you have an idea how to use RandomPoint[Sphere[], n] for regular self-intersecting random walk.

Graphics3D[
    {{Opacity[.5],s},
    {PointSize[.005],Opacity[.5],Point[pts]},
    {Opacity[.7],Line[walk]}}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.