5
$\begingroup$

I would like to plot a sphere with multi-colored regions separated by solutions of a differential equation (no explicit expression is available for the separating curves).

The result I am after would look similar to the one produced by

 Show[ParametricPlot3D[{Sin[t] Cos[p], Sin[t] Sin[p], Cos[t]}, {t, 0, π}, 
       {p, 0, 2 π}, ColorFunction -> Function[{x, y, z},        
        If[-0.3 < z < 0.3, {Red, Opacity[0.8]}, {Green, Opacity[0.8]}]], 
      ColorFunctionScaling -> False, PlotPoints -> 300, Mesh -> 0], 
     ParametricPlot3D[{Sqrt[0.91] Sin[p], Sqrt[0.91] Cos[p], 0.3}, {p, 0, 
       2 π}, PlotStyle -> Directive[Thick, Black], PlotPoints -> 100],
      ParametricPlot3D[{Sqrt[0.91] Sin[p], Sqrt[0.91] Cos[p], -0.3}, {p, 
       0, 2 π}, PlotStyle -> Directive[Thick, Black], 
      PlotPoints -> 100]]

Colored regions

except the lines separating the green and red regions in my case need to be the solutions of a system of differential equations:

a = 1;
b = 5;
c = 20;
T = 10;
p1 = {-0.8, 0., -0.6};
p2 = {0.8, 0., 0.6};
sol1 = NDSolve[{x'[t] == y[t] (b x[t] - c z[t]), 
    y'[t] == -b x[t]^2 + c x[t] z[t] - z[t] (a + b z[t]), 
    z'[t] == y[t] (a + b z[t]), {x[0], y[0], z[0]} == p1}, {x[t], 
    y[t], z[t]}, {t, 0, T}];
sol2 = NDSolve[{x'[t] == y[t] (b x[t] - c z[t]), 
    y'[t] == -b x[t]^2 + c x[t] z[t] - z[t] (a + b z[t]), 
    z'[t] == y[t] (a + b z[t]), {x[0], y[0], z[0]} == p2}, {x[t], 
    y[t], z[t]}, {t, 0, T}];
plot2 = Show[
  ParametricPlot3D[{Sin[t] Cos[p], Sin[t] Sin[p], Cos[t]}, {t, 
    0, π}, {p, 0, 2 π}, PlotStyle -> Directive[Opacity[0.8]], 
   PlotPoints -> 100, Mesh -> 0], 
  ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol1], {t, 0, T}, 
   PlotStyle -> Directive[Black, Thick]], 
  ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol2], {t, 0, T}, 
   PlotStyle -> Directive[Black, Thick]]]

How do I use the InterpolatingFunction trajectories sol1 and sol2 from NDSolve to produce multi-colored regions on a sphere? Regions to color

$\endgroup$
1
  • $\begingroup$ May we constrain the problem to two non-intersecting closed curves on the surface of the sphere? $\endgroup$
    – m_goldberg
    Sep 10, 2015 at 1:27

1 Answer 1

9
$\begingroup$

Fortunately, your solutions look very friendly, because I can easily think of situations where the creation of such a plot would be harder. Let me try to illustrate how you can solve one part of your problem and leave the rest for yourself. In the end it should look like this

Mathematica graphics

What we have are the solutions to your differential equation in Cartesian coordinates. I added a First to the call of NDSolve, the rest is unchanged:

a = 1;
b = 5;
c = 20;
T = 10;
p1 = {-0.8, 0., -0.6};
p2 = {0.8, 0., 0.6};
sol1 = First@
   NDSolve[{x'[t] == y[t] (b x[t] - c z[t]), 
     y'[t] == -b x[t]^2 + c x[t] z[t] - z[t] (a + b z[t]), 
     z'[t] == y[t] (a + b z[t]), {x[0], y[0], z[0]} == p1}, {x, y, 
     z}, {t, 0, T}];
sol2 = First@
   NDSolve[{x'[t] == y[t] (b x[t] - c z[t]), 
     y'[t] == -b x[t]^2 + c x[t] z[t] - z[t] (a + b z[t]), 
     z'[t] == y[t] (a + b z[t]), {x[0], y[0], z[0]} == p2}, {x, y, 
     z}, {t, 0, T}];

I'm not sure whether you have ever thought about how ParametricPlot3D draws a sphere, so let us start with this:

    ParametricPlot3D[{Sin[θ] Cos[ϕ], Sin[θ] Sin[ϕ], Cos[θ]}, 
      {θ, 0, π}, {ϕ, -Pi, Pi}]

What basically happens is that {Cos[ϕ],Sin[ϕ]} draws a circle in the x-y-plane. This is why ϕ runs through the whole 2Pi interval. To make a complete sphere, this ring is now moved up- and downwards in z-direction which is why there is a Cos[θ] in the z-coordinate. Additionally, the ring needs to get smaller because otherwise we would end in a simple cylinder. The scaling is done, by multiplying the ring with Sin[θ].

The important part is to understand that we have two directions. We move along the circle-ring with ϕ and we move up- and downwards with θ. What if we use instead of a circle, your ring like structure in this approach?

First problem is that we don't have your solution in spherical angles, but in Cartesian coordinates. So let's tackle this first. Look at this

CoordinateTransform[
 "Cartesian" -> "Spherical", {x[t], y[t], z[t]}]
(* {Sqrt[x[t]^2 + y[t]^2 + z[t]^2], 
 ArcTan[z[t], Sqrt[x[t]^2 + y[t]^2]], ArcTan[x[t], y[t]]} *)

Using this transformation, we can create functions that give use the corresponding angles (and r) when we input a specific t

r[t_?NumericQ] = Sqrt[x[t]^2 + y[t]^2 + z[t]^2] /. sol1;
theta[t_?NumericQ] = ArcTan[z[t], Sqrt[x[t]^2 + y[t]^2]] /. sol1;
phi[t_?NumericQ] = ArcTan[x[t], y[t]] /. sol1;
Plot[{r[t], theta[t], phi[t]}, {t, 0, 1}, 
 PlotLegends -> {r, theta, phi}, PlotRange -> All]

Mathematica graphics

The first thing that jumps directly into our eye is that your function is periodic at a t-value of about 0.5 (and not your T of 10). When we will draw your ring, we don't want to draw them over and over again, because this ends usually in visual artifacts. Especially when we deal with surfaces. You can find the correct end time by

tEnd = t /. FindRoot[phi[0] - phi[t] /. sol1, {t, .4, .6}]
(* 0.496357 *)

Now we are all set. What we will do now is to not iterate the angles directly, but we iterate t and compute the angles with the functions phi and theta we have created. We won't need r btw. Therefore, the first step to solve your problem is

Show[
 ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol1], {t, 0, tEnd}, 
   PlotStyle -> Directive[Gray]] /. Line[pts__] :> Tube[pts, 0.01],

 ParametricPlot3D[{Sin[theta2] Cos[phi[t]], Sin[theta2] Sin[phi[t]], 
   Cos[theta2]},
  {t, 0, tEnd}, {theta2, theta[t], Pi}, 
  PlotStyle -> Directive[Opacity[0.5, Green]], Mesh -> None, 
  MaxRecursion -> 3],

 ParametricPlot3D[{Sin[theta2] Cos[phi[t]], Sin[theta2] Sin[phi[t]], 
   Cos[theta2]},
  {t, 0, tEnd}, {theta2, 0, theta[t]}, 
  PlotStyle -> Directive[Opacity[0.8, Red]], Mesh -> None, 
  MaxRecursion -> 3],
 PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}
 ]

Mathematica graphics

The first parametric plot creates the ring. The second one draws the small cup by iterating theta2 from the current angle that points on the ring theta[t] to the very bottom at 0. The other cup is created by exactly the same call only that theta2 iterates over the other half.

I have left the big part for yourself to not spoil the fun: drawing the middle part. The other smaller cup should be similar to what I explained here.

Import["http://goo.gl/NaH6rM"]["https://i.stack.imgur.com/fouCb.png"]

$\endgroup$
1
  • $\begingroup$ Thank you, halirutan! I was able to achieve what I wanted. $\endgroup$ Sep 10, 2015 at 22:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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