Hello
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
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]
plt
I am getting this plot-
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
g(x,y,z)=-x+y+z
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)]
end
#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):
n=200
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",
))
end
display(fig)
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:
NB:
The result is not as expected, as per @sijoâs post below. Needs more work.
Rome ne sâest pas faite en un jourâŠ
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)
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âŠ
@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):
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)
end
end
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:
fvals
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
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: https://github.com/empet/Datasets/blob/master/Images/sol-linineq1.png
