# How to plot multi-colored regions separated by InterpolatingFunction on sphere

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]]


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?

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

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

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]


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}}
]


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"]

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