# why ContourPlot3D does not give plot?

I wanted to plot the intersections using ContourPlot3D in order to observe how do solutions look like, however, the following code does not work:

ContourPlot3D[{Sin[a - b] + Sin[a] + Sin[a - d] == 0 &&
Sin[b - a] + Sin[b] + Sin[b - d] == 0 &&
Sin[-a] + Sin[-b] - Sin[-d] == 0 &&
Sin[d - a] + Sin[d - b] == 0}, {a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0,
2 Pi}, PlotPoints -> 50]

Could anyone help? Thanks in advance.

• Replace && with commas to get this. Increase the points only if you think it is worth it.
– Syed
Feb 7, 2022 at 19:59
• Alternatively, you can plot individually and then call them collectively, like so: plot1 = ContourPlot3D[{Sin[a - b] + Sin[a] + Sin[a - d] == 0}, {a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi}]; plot2 = ContourPlot3D[{Sin[b - a] + Sin[b] + Sin[b - d] == 0}, {a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi}]; plot3 = ContourPlot3D[{Sin[-a] + Sin[-b] - Sin[-d] == 0}, {a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi}]; plot4 = ContourPlot3D[{Sin[d - a] + Sin[d - b] == 0}, {a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi}]; and to make them appear all together Show[plot1, plot2, plot3, plot4]
– user49048
Feb 7, 2022 at 20:03
• In 3D, 1 equation will give a surface, 2 equations give a 1 dimensional region that can be empty. 3 equations may give a 0 dimensional region, that may also be empty. However, in general, 4 equations do not have a solution. Feb 7, 2022 at 20:05
• @Seyd. Replacing && by commas will plot all the surfaces belonging to the different single equations and not their intersection. Feb 7, 2022 at 20:09
• @DanielHuber Indeed, it is all the surfaces instead of intersection. Could I ask: why in general 4 equations do not have a solution? Thank you!
– M.K
Feb 7, 2022 at 20:13

We can simplify the original equations by GroebnerBasis method.

Clear["Global*"];
trigs = {Sin[-a] + Sin[-b] - Sin[-d],
Sin[a - b] + Sin[a] + Sin[a - d], Sin[b - a] + Sin[b] + Sin[b - d],
Sin[d - a] + Sin[d - b]};
polys = TrigExpand /@ trigs /. {Sin[a] -> sa, Sin[b] -> sb,
Sin[d] -> sd, Cos[a] -> ca, Cos[b] -> cb, Cos[d] -> cd};
allfuns =
Join[polys, {sa^2 + ca^2 - 1, sb^2 + cb^2 - 1, sd^2 + cd^2 - 1}];
gb = GroebnerBasis[allfuns, {sa, sb, sd}, {ca, cb, cd}];
sol = Reduce[gb == 0] /. {sa -> Sin[a], sb -> Sin[b], sd -> Sin[d]};
oldeqns = And @@ Thread[trigs == 0];
neweqns = FullSimplify[oldeqns && sol]
Simplify[neweqns, oldeqns];
Simplify[oldeqns, neweqns];

Sin[a - b] == Sin[b] + Sin[b - d] && Sin[a - d] + Sin[b - d] == 0 && Sin[d] == 0 && Sin[a] + Sin[b] == 0

result = Reduce[
neweqns && 0 <= a <= 2 π && 0 <= b <= 2 π &&
0 <= d <= 2 π, {a, b, d}, Reals] // FullSimplify

reg = ImplicitRegion[neweqns, {a, b, d}];
instances =
FindInstance[{a, b, d} ∈ reg && 0 <= a <= 2 π &&
0 <= b <= 2 π && 0 <= d <= 2 π, {a, b, d}, 1000];
AllTrue[oldeqns /. instances, TrueQ]
Graphics3D[{Red, Point[{a, b, d} /. instances]}]

• Thank you! What is 'Simplify[neweqns, oldeqns]' and 'Simplify[oldeqns, neweqns]' trying do?
– M.K
Feb 9, 2022 at 14:53
• @HJ_dynamics It means that neweqns is equivalent to oldeqns. Feb 9, 2022 at 15:16
Clear["Global*"]

eqns =
{Sin[a - b] + Sin[a] + Sin[a - d] == 0,
Sin[b - a] + Sin[b] + Sin[b - d] == 0,
Sin[-a] + Sin[-b] - Sin[-d] == 0,
Sin[d - a] + Sin[d - b] == 0};

Solving within the PlotRange

sol = Solve[Join[eqns, Thread[0 <= {a, b, d} <= 2 Pi]],
{a, b, d},
Method -> Reduce,
MaxExtraConditions -> 1] // Simplify

(* Solve::fexp: Warning: Solve used FunctionExpand to transform the system. Since FunctionExpand transformation rules are only generically correct, the solution set might have been altered.

Solve::svars: Equations may not give solutions for all "solve" variables.

{{b -> ConditionalExpression[a - π,
a ∈ Reals && π <= a <= 2 π],
d -> ConditionalExpression[π,
a ∈ Reals && π <= a <= 2 π]}, {b ->
ConditionalExpression[a + π, a ∈ Reals && 0 <= a <= π],
d -> ConditionalExpression[π,
a ∈ Reals && 0 <= a <= π]}, {a -> 0, b -> 0,
d -> 0}, {a -> 0, b -> 0, d -> π}, {a -> 0, b -> 0,
d -> 2 π}, {a -> 0, b -> π, d -> 0}, {a -> 0, b -> π,
d -> 2 π}, {a -> 0, b -> 2 π, d -> 0}, {a -> 0, b -> 2 π,
d -> π}, {a -> 0, b -> 2 π, d -> 2 π}, {a -> π, b -> 0,
d -> 0}, {a -> π, b -> 0, d -> 2 π}, {a -> π, b -> π,
d -> 0}, {a -> π, b -> π, d -> π}, {a -> π, b -> π,
d -> 2 π}, {a -> π, b -> 2 π, d -> 0}, {a -> π,
b -> 2 π, d -> 2 π}, {a -> 2 π, b -> 0, d -> 0}, {a -> 2 π,
b -> 0, d -> π}, {a -> 2 π, b -> 0, d -> 2 π}, {a -> 2 π,
b -> π, d -> 0}, {a -> 2 π, b -> π,
d -> 2 π}, {a -> 2 π, b -> 2 π, d -> 0}, {a -> 2 π,
b -> 2 π, d -> π}, {a -> 2 π, b -> 2 π,
d -> 2 π}, {a -> π/2, b -> (3 π)/2,
d -> π}, {a -> (3 π)/2, b -> π/2, d -> π}} *)

Verifying the solutions,

(And @@ eqns) /. sol

(* {ConditionalExpression[True, a ∈ Reals && π <= a <= 2 π],
ConditionalExpression[True,
a ∈ Reals &&
0 <= a <= π], True, True, True, True, True, True, True, True, True, \
True, True, True, True, True, True, True, True, True, True, True, True, True, \
True, True, True} *)

Plotting,

Show[
ContourPlot3D[Evaluate@eqns,
{a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi},
AxesLabel -> (Style[#, 14] & /@ {a, b, d}),
PlotLegends -> "Expressions"],
Graphics3D[{Red, AbsolutePointSize[8],
Point[{a, b, d} /. sol[[3 ;;]]]}]]

Show[
ContourPlot3D[Evaluate@eqns,
{a, 0, 2 Pi}, {b, 0, 2 Pi}, {d, 0, 2 Pi},
ContourStyle -> Opacity[0.05],
Mesh -> None,
AxesLabel -> (Style[#, 14] & /@ {a, b, d})],
Graphics3D[{Red, AbsolutePointSize[8],
Point[{a, b, d} /. sol[[3 ;;]]]}]]