I am trying to solve a stochastic equation corresponding to Vector resonant relaxation (the way orbital planes of stars diffuse near the galactic center, see below).

enter image description here


Within the context of some approximation, the orientation vector of the plane of the orbit $\hat{\mathbf{L}}$ obeys

enter image description here

where $\mathbf{M}$ is a stochastic matrix of the form

enter image description here

which I have implemented in Mathematica as (assuming $\Gamma_t=1$)

Format[em2] = Subscript[e, -2]; Format[em1] = Subscript[e, -1];
Format[e2] = Subscript[e, 2];
Format[e1] = Subscript[e, 1];
Format[e0] = Subscript[e, 0]; Format[u1] = Subscript[u, 1]; 
Format[u2] = Subscript[u, 2]; Format[u3] = Subscript[u, 3];

M = {{DifferentialD[e2][t], 
   DifferentialD[em2][t], -DifferentialD[e1][t]}, {DifferentialD[em2][
    t], -DifferentialD[e2][t], -DifferentialD[em1][
     t]}, {-DifferentialD[e1][t], -DifferentialD[em1][
     t], Sqrt[3] DifferentialD[e0][t]}}

Then I can define $\hat{ \mathbf{L}}$ as

L = #[t] & /@ {u1, u2, u3};
eqn = \[DifferentialD]#[t] & /@ {u1, u2, u3} == Cross[L, M.L] // Apart // Thread;
eqn // TableForm

enter image description here

So that the ItoProcess obeys

proc = ItoProcess[
  eqn, {u1[t], u2[t], u3[t]}, {{u1, u2, u3}, {1, 1/2, 1/3}}, t,
  {e2 \[Distributed] WienerProcess[], 
   em2 \[Distributed] WienerProcess[], 
   e1 \[Distributed] WienerProcess[], 
   em1 \[Distributed] WienerProcess[], 
   e0 \[Distributed] WienerProcess[]}]

enter image description here

I can then compute a set of paths as

path = RandomFunction[proc, {0., 1., 0.01}, Method -> "StochasticRungeKutta"]

enter image description here


Why is it not returning random paths?

Let me apologize in advance for what is likely to be a stupid mistake on my part.

Astrophysical context

In the vicinity of a supermassive black hole, stars move on nearly Keplerian orbits. Yet, because of the enclosed stellar mass and relativistic corrections, the potential slightly deviates from the Keplerian one, which causes the stellar orbits to precess. Similarly, as a result of the finite number of stars, the mutual gravitational torques between pairs of stars also drive a rapid reshuffling of the stars' orbital orientations. Overall, the combination of these two effects leads to a stochastic evolution of stellar orbital angular momentum vectors, through a process named 'resonant relaxation'.


The typographic error is in the definition of the matrix

M = {{DifferentialD[e2[t]], 
   DifferentialD[em2[t]], -DifferentialD[e1[t]]}, {DifferentialD[
   -DifferentialD[e2[t]], -DifferentialD[em1[t]]}, {-DifferentialD[
     e1[t]], -DifferentialD[em1[t]], Sqrt[3] DifferentialD[e0[t]]}}

i.e. it should involve DifferentialD[e2[t]] or \[DifferentialD]e2[t] not DifferentialD[e2][t].

enter image description here

Since for this problem the noise needs to be correlated (see this question), it is in fact simpler to rely on NDSolve together with GaussianRandomField

nn = 1024*2; tmax = 2*25; Amp = 5;
noise = Interpolation[#, Method -> "Spline"][t] & /@ 
   Table[Transpose@{tmax Range[0, nn - 1]/(nn - 1) // N,
      GaussianRandomField[nn, 1, Function[k, k^-0.5]] Amp}, {5}];
eqn2 = eqn /. {\[DifferentialD]e0[t] -> 
      noise[[1]], \[DifferentialD]e1[t] -> 
      noise[[2]], \[DifferentialD]em1[t] -> 
      noise[[3]], \[DifferentialD]e2[t] -> 
      noise[[4]], \[DifferentialD]em2[t] -> 
      noise[[5]]} /.
   {\[DifferentialD]u1[t] -> 
     u1'[t], \[DifferentialD]u2[t] -> 
     u2'[t], \[DifferentialD]u3[t] -> u3'[t]};
eqn2 = Join[eqn2, {u1[0] == 1, u2[0] == 0, u3[0] == 0}];
sol[t_] = NDSolveValue[eqn2, {u1[t], u2[t], u3[t]}, {t, 0, tmax}]


Show[ParametricPlot3D[ sol[t], {t, 0, tmax}, BoxRatios -> {1, 1, 1}], 
 Graphics3D[{Opacity[0.2], Sphere[]}], Boxed -> False, Axes -> False]

enter image description here

One can represent explicitly the orbital planes via this function

circle3D[centre_: {0, 0, 0}, radius_: 1, normal_: {0, 0, 1}, 
  angle_: {0, 2 Pi}] := 
   Map[RotationTransform[{{0, 0, 1}, normal}, centre], #] &, 
   Map[Append[#, Last@centre] &, #] &, 
   Append[DeleteDuplicates[Most@#], Last@#] &, Level[#, {-2}] &, 
   MeshPrimitives[#, 1] &, DiscretizeRegion, If][
  First@Differences@angle >= 2 Pi, Circle[Most@centre, radius], 
  Circle[Most@centre, radius, angle]]

from that answer. Then

   Table[circle3D[{0, 0, 0}, i, sol[i]], {i, 0, 1, 0.05}]}, 
  Boxed -> False], Boxed -> False, Axes -> False]

enter image description here


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.