How to plot a 3D plot of 3 inequalities?


I have 3 inequalities and each have 3 variables. Could you please explain how can I draw a combined one 3D plot with shaded area satisfying all the 3 inequalities? I have inequalities of the form as

f1(x,y,z) = a1*x + b1*y+ c1*z <= 0
f2(x,y,z) = a2*x + b2*y+ c2*z <= 0
f3(x,y,z) = a3*x + b3*y+ c3*z <= 0

Thank you

A simple solution, but in the form of a 3d scatter plot is given here.

Dear @rafael.guerra,
Thank you for providing a solution. But I am trying to plot using “IntervalArithmetic.jl” package.

I tried to plot using the code shown below-

f1(p, q, r) = p^2 + q^2 + r^2 - 4 #sphere
	f2(p, q, r) = p^2 + q^2 + r^2 - 16  #sphere
	# f3(p, q, r) = 2*p + 3*q
	zrange = -5:.1:5
	plt = plot(palette = palette(cgrad([:red,:green,:blue],length(zrange))),camera=(45,45),xlims=(-7,7), ylims=(-7,7), xlabel="x",ylabel="y",zlabel="z")
	[plot!(plt,((p, q)->f1(p, q, r)) ≧ 0, zpos=r) for r=zrange]
	[plot!(plt,((p, q)->f2(p, q, r)) ≦ 0, zpos=r) for r=zrange]
	# [plot!(plt,((p, q)->f3(p, q, r)) ⩵ 0, zpos=r) for r=zrange]

I am getting this plot-
Image 02-02-23 at 4.30 PM

But, don’t know why these additional square boxes and and lines are coming. I just want spherical region in the plot. Could you please help how can I remove everything other than the spherical region.

Thank you

Your region of interest is delimited by three planes, passing all through the origin. Hence you should plot the three planes as parameterized surfaces. I don’t know how this could be plotted with Plots, but I explain with a PlotlyJS example.
Suppose that the three planes are defined by the linear functions:

f(x,y,z) =2x-y+3z
h(x,y,z) =x-3y-5z
using LinearAlgebra, PlotlyJS
import StaticArrays: SMatrix
function planeparam(normal::SMatrix{1, 3, T}; p=[0,0,0]) where T<:AbstractFloat
    ~iszero(normal) || error("normal vector is null vect")
    M = nullspace(normal)
    a, b = eachcol(M)
    #parameterization of the plane π(p; normal=no) defined by the point p and the given normal
    (u,v)->[a[k]*u+b[k]*v+p[k] for k in 1:3]  # theoretically: (u, v)->[X(u, v), Y(u,v), Z(u,v)]

#Get the parameterization of each plane
F = planeparam(SMatrix{1,3}([2.0 -1.0 3]); )
G = planeparam(SMatrix{1,3}([-1.0  1  1.]); )
H = planeparam(SMatrix{1,3}([1.0  -3  -5.]); )
#define a meshgrid over a rectangular region in the parametric plane (u,v):
ul = vl= range(-3, 3, n)
u = ones(n) * ul'
v = vl *ones(n)'

fig= Plot(Layout(width=450, height=450, font_size=11,
                 coloraxis = attr(colorscale=colors.algae, reversescale=true,  
                                  showscale=false ),
                scene_camera_eye=attr(x=1.5, y=1.5, z=1)))
#add the planes as parameterized surfaces to the figure fig
for K in [F, G, H]
    X, Y, Z = eachrow(hcat(K.(u,v)...))
    addtraces!(fig, surface( x=reshape(X, n, n), y=reshape(Y, n, n), 
               z=reshape(Z, n, n), opacity=0.75, coloraxis="coloraxis", 


I selected three random planar surfaces, hence the plot is not quite convincing, but it can suggest an idea on how to proceed. Here is missing the condition <=0. For linear functions it is very simple to get the region whree both three are negative, but it is complicated to shade that volume, because it is covered by the surfaces.
Note that depending on your linear functions, the volume where all three are negative can be bounded or unbounded.


This problem seems to correspond exactly to Polyhedra.jl’s H-representation of halfspaces.

We can use Polyhedra.jl and Makie.jl to plot it:

The result is not as expected, as per @sijo’s post below. Needs more work.
Rome ne s’est pas faite en un jour…

Polyhedra.jl and Makie.jl code
using Polyhedra
using Makie, GLMakie

const a1, b1, c1 = 0.,  1., 1.
const a2, b2, c2 = 1.,  1., 1.
const a3, b3, c3 = 1., -1., 1.

f1(x,y,z) = a1*x + b1*y + c1*z
f2(x,y,z) = a2*x + b2*y + c2*z
f3(x,y,z) = a3*x + b3*y + c3*z

h = HalfSpace([a1, b1, c1], 0) ∩ HalfSpace([a2, b2, c2], 0) ∩ HalfSpace([a3, b3, c3], 0)
p = polyhedron(h)
m = Polyhedra.Mesh(p)

Makie.mesh(m, color=:blue)
1 Like

It’s a good idea to use Polyhedra.jl when the solution is a bounded volume, but if it unbounded it throws the error:

 hred: Redundancy Polyhedra.UNKNOWN_REDUNDANCY
  vrep: Nothing nothing
  vred: Redundancy Polyhedra.UNKNOWN_REDUNDANCY
  solver: Nothing nothing

But maybe for OP’s case it is bounded.

I don’t know what Polyhedra.jl is doing in @rafael.guerra’s example but you can’t have a (nonzero) bounded volume delimited by three planes…

1 Like

@sijo, good observation and I do not know the answer either.

Using the brute force scatter3d method, as per the Python post linked previously, we get an unbounded volume (computed over limited range here):

Plots.jl plotly() code
const a1, b1, c1 = 0.,  1., 1.
const a2, b2, c2 = 1.,  1., 1.
const a3, b3, c3 = 1., -1., 1.

f1(x,y,z) = a1*x + b1*y + c1*z
f2(x,y,z) = a2*x + b2*y + c2*z
f3(x,y,z) = a3*x + b3*y + c3*z

sup = -10:0.2:10
X = Float32[]
Y = Float32[]
Z = Float32[]
for x in sup, y in sup, z in sup
    if f1(x,y,z) ≤ 0 && f2(x,y,z) ≤ 0 && f3(x,y,z) ≤ 0
        push!(X, x)
        push!(Y, y)
        push!(Z, z)

using Plots; plotly(size=(600,600))
scatter3d(X,Y,Z, ms=1, msw=0.1, xlabel="x", ylabel="y", zlabel="z", c=:blues)

Dear @rafael.guerra,
As I don’t know whether the actual volume will be bounded or unbounded, this type of solution (plot) seems useful to me. Could you please share julia code and julia packages used for this?
Thank you

Included it below the plot.

Update: Such a system of linear inequalities can be solved with MarchingCubes.jl as follows:

  1. define a vector valued function f(x,y, z) =[f1(x,y,z), f2(x,y,z), f3(x,y,z)],
    and a cube with center at origin, \text{Cube}= [-a, a]^3
  2. a grid on this cube and a function that gets 0 at the boundary of the region of interest, i.e. the volume V=\{(x,y,z)\:|\: f1(x,y, z)<=0, f2(x,y, z)<=0 , f3(x,y, z)<=0 \}. Evaluate this function at the grid points, to get fvals
  3. Call the marching cube method to retrieve a triangulation of the volume boundary and plot the corresponding mesh.
using MarchingCubes
f(x,y,z) = [y + z, x+y+z, x-y+z]
#f(x,y, z) = [2*x-4.5*y+3.2*z, -3x+2.1*y+1.0*z, x-3*y-5*z] 
const Fl=Float64
#Define a cube  of side 14
xl = yl = zl = collect(Float64,range(-7.0, 7.0, length=150))
x = [xx for zz in zl, yy in yl, xx in xl]
y = [yy for zz in zl, yy in yl, xx in xl]
z = [zz for zz in zl, yy in yl, xx in xl];
fvals = [maximum([f(xx, yy, zz)[k] for k=1:3]) for zz in zl, yy in yl,  xx in xl]
#instantiate the structure MC (Marching Cube)
mc = MC(fvals, Int; x=xl, y=yl, z=zl)
m = march(mc, 0.0) #0.0 is the isovalue
  1. Transform mc.vertices and mc.triangles to data types necessary to plot a 3d mesh in your favorite graphics library.

This is the plot I created with PlotlyJS.mesh3d: