# How can I find the intersection of two sets and how can I sample elements from this intersection?

I've solved a non-linear system of equations and the numerical solution I obtained was the following:

  {{x ->
ConditionalExpression[0.125 (-5. - 8. y + 48. z),
(-24.0367 <= y <= 0.373823
&&
Root[1.4184*10^9 - 1.14866*10^9 y + 3.59484*10^8 y^2 -
1.83744*10^8 y^3 +
5.568*10^7 y^4 + (-1.65159*10^10 - 5.23081*10^9 y -
6.91386*10^9 y^2 - 9.78854*10^8 y^3 -
6.5257*10^8 y^4) #1 + (-9.28089*10^9 - 9.79188*10^8 y -
1.24545*10^9 y^2 + 1.60358*10^8 y^3 +
1.06906*10^8 y^4) #1^2 + (-2.31117*10^8 +
3.12699*10^8 y + 4.16932*10^8 y^2) #1^3 +
4.26553*10^8 #1^4 &, 1]
<= z <=
Root[1.77292*10^11 - 1.43567*10^11 y + 4.49381*10^10 y^2 -
2.29657*10^10 y^3 +
6.96*10^9 y^4 + (-2.06444*10^12 - 6.53839*10^11 y -
8.64219*10^11 y^2 - 1.22355*10^11 y^3 -
8.15703*10^10 y^4) #1 + (-1.1601*10^12 -
1.22397*10^11 y - 1.55679*10^11 y^2 +
2.00448*10^10 y^3 +
1.33632*10^10 y^4) #1^2 + (-2.88882*10^10 +
3.90874*10^10 y + 5.21165*10^10 y^2) #1^3 +
5.33192*10^10 #1^4 &, 2])
||
(Root[1.77292*10^11 - 1.43567*10^11 y + 4.49381*10^10 y^2 -
2.29657*10^10 y^3 +
6.96*10^9 y^4 + (-2.06444*10^12 - 6.53839*10^11 y -
8.64219*10^11 y^2 - 1.22355*10^11 y^3 -
8.15703*10^10 y^4) #1 + (-1.1601*10^12 -
1.22397*10^11 y - 1.55679*10^11 y^2 +
2.00448*10^10 y^3 +
1.33632*10^10 y^4) #1^2 + (-2.88882*10^10 +
3.90874*10^10 y + 5.21165*10^10 y^2) #1^3 +
5.33192*10^10 #1^4 &, 1]
<= z <=
Root[1.77292*10^11 - 1.43567*10^11 y + 4.49381*10^10 y^2 -
2.29657*10^10 y^3 +
6.96*10^9 y^4 + (-2.06444*10^12 - 6.53839*10^11 y -
8.64219*10^11 y^2 - 1.22355*10^11 y^3 -
8.15703*10^10 y^4) #1 + (-1.1601*10^12 -
1.22397*10^11 y - 1.55679*10^11 y^2 +
2.00448*10^10 y^3 +
1.33632*10^10 y^4) #1^2 + (-2.88882*10^10 +
3.90874*10^10 y + 5.21165*10^10 y^2) #1^3 +
5.33192*10^10 #1^4 &, 2]
&&
y > 0.373823)
||
(Root[1.77292*10^11 - 1.43567*10^11 y + 4.49381*10^10 y^2 -
2.29657*10^10 y^3 +
6.96*10^9 y^4 + (-2.06444*10^12 - 6.53839*10^11 y -
8.64219*10^11 y^2 - 1.22355*10^11 y^3 -
8.15703*10^10 y^4) #1 + (-1.1601*10^12 -
1.22397*10^11 y - 1.55679*10^11 y^2 +
2.00448*10^10 y^3 +
1.33632*10^10 y^4) #1^2 + (-2.88882*10^10 +
3.90874*10^10 y + 5.21165*10^10 y^2) #1^3 +
5.33192*10^10 #1^4 &, 1] <= z <=
Root[1.77292*10^11 - 1.43567*10^11 y + 4.49381*10^10 y^2 -
2.29657*10^10 y^3 +
6.96*10^9 y^4 + (-2.06444*10^12 - 6.53839*10^11 y -
8.64219*10^11 y^2 - 1.22355*10^11 y^3 -
8.15703*10^10 y^4) #1 + (-1.1601*10^12 -
1.22397*10^11 y - 1.55679*10^11 y^2 +
2.00448*10^10 y^3 +
1.33632*10^10 y^4) #1^2 + (-2.88882*10^10 +
3.90874*10^10 y + 5.21165*10^10 y^2) #1^3 +
5.33192*10^10 #1^4 &, 2]
&&
y < -24.0367)]}}


How can I find the intersection of two sets and how can I sample elements from this intersection?

To be honest, I can't even figure out what elements belong to the set I obtained as the solution to my system of equations. But that's OK, because all I'm interested in is finding its intersection with the set

{x^2 + y^2 + z^2 <= 1} - {x^2 + y^2 + z^2 <= 0.001}.


Edit 3: The set I mean above is

{0.001 <= x^2 + y^2 + z^2 <= 1}

I'd also like to sample elements from this intersection in case it is non-empty. Is it possible to do such things with Mathematica?

Edit 1: Here follows the code that gave me that solution:

Edit 2: I changed the function k

k = (((6 - x)^2 - (y)^2)*z + 0.5*(6^2 - (6 - x)^2)*10 +
0.5*10*y^2)/(2*z*(6 - x - y) + x*10 + y*10)

Ixx = 1/3*(10*((6 - (k))^3 - (6 - x - (k))^3) +
2*z*((6 - x - (k))^3 - (-(k) + y)^3) +
10*((-(k) + y)^3 - (-(k))^3))

Q1 = -(k)*(y*10 + 2*((6 - x) - (k))*z) + 0.5*10*(6^2 - (6 - x)^2) +
z*((6 - x)^2 - (k)^2)

Q2 = (k)*(y*10 + 2*((k) - y)*z) + 0.5*10*y^2 - z*((k)^2 - (y)^2)
NSolve[{290*10*Ixx >= 2.4*Q1, 580*Ixx >= 24*6/2,
60*(2*6*(z - x - y) + x*10 + y*10)*8*10 + 2*500*12 ==
60*(10*3)*1*10}, {x, y, z}, Reals]

• show you initial code of before solution! to make our life easier ! Apr 4, 2019 at 23:35
• @Alrubaie sure thing! Apr 5, 2019 at 2:24

We can avoid the conditional expression altogether with the following approach. First, your expressions

ClearAll["Global*"]

k = (6^2*z + 0.5*(6 - x)*10 + 0.5*10*y^2)/(2*6*(z - x - y) + x*10 +
y*10);

Ixx = 1/3*(10*((6 - (k))^3 - (6 - x - (k))^3) +
2*z*((6 - x - (k))^3 - (-(k) + y)^3) +
10*((-(k) + y)^3 - (-(k))^3));

Q1 = -(k)*(y*10 + 2*((6 - x) - (k))*z) + 0.5*10*(6^2 - (6 - x)^2) +
z*((6 - x)^2 - (k)^2);

Q2 = (k)*(y*10 + 2*((k) - y)*z) + 0.5*10*y^2 - z*((k)^2 - (y)^2);


Now, instead of using NSolve we can write the two inequalities and the one equation as

c1 = 290*10*Ixx >= 2.4*Q1;
c2 = 580*Ixx >= 24*6/2;

eqn = 60*(2*6*(z - x - y) + x*10 + y*10)*8*10 + 2*500*12 ==
60*(10*3)*1*10;


We can obtain a value, zeta, for z using Solve, as

ζ = z /. First@Solve[eqn, z]


We can now define a region in the (x,y)-plane that satisfies our two inequality constraints, like this

r1 = ImplicitRegion[ c1 && c2 /. z->ζ,
{{x,-10,50},{y,-10,10}}];
Plot3D[ζ, {x, y} ∈ r1]


The result above is a visualization of the region of $$z=\zeta$$ plane that satisfies the inequalities and the equation that had been in the NSolve expression.

We have two more constraints that we can use. If we only want to consider solutions between the two spheres, we can write these expressions

c3 = {x, y, z} ∈ Ball[{0, 0, 0}, 1];
c4 = {x, y, z} ∈ Ball[{0, 0, 0}, 1/Sqrt[1000]];

r2 = ImplicitRegion[c1 && c2 && c3 /. z -> ζ,
{{x, -10, 50}, {y, -10, 10}}];
Plot3D[{ζ, Sqrt[1/1000 - x^2 - y^2]}, {x, y} ∈ r2]


That gives us a visualization of the region in the solution between the two spheres and of the smaller sphere.

Note that the second region, $$r_2$$, does not explicitly include && Not[c4] in its constraints. We can check the minimum $$z$$ value for the $$(x,y)$$ points in $$r_2$$ with

Minimize[{ζ, {x, y} ∈ r2}, {x, y}]
(*  {0.0427734, {x -> -0.866208, y -> 0.497849}}  *)


This tell us that the region $$z=\zeta$$ does not intersect the ball $$c_4$$ for $$(x,y)\in r_2$$.

Finally, if we want to sample the points in $$r_2$$, we can use this:

nSample = 10;
{x, y, ζ} /. FindInstance[{x, y} ∈ r2, {x, y}, nSample]

• Thank you very much for your answer, really appreciate it! Would you know how could I restrict the solution set to values of x,y,z>0? Apr 5, 2019 at 13:30
• We can restrict the sampling region to positive values of x, y, z by modifying the intervals in the ImplicitRegion` statement that defines the region. Apr 5, 2019 at 18:47
• I just realized something, when I said I was only interested in the solutions satisfying {x^2 + y^2 + z^2 <= 1} - {x^2 + y^2 + z^2 <= 0.001}, I meant to say those that are elements of the set {0.001 <= x^2 + y^2 + z^2 <= 1}. I tried correcting that on your code by defining c4 as: c4 = {x, y, z} [NotElement] Ball[{0, 0, 0}, 0.001]; And redefining c2 as: r2 = ImplicitRegion[ c1 && c2 && c3 && c4 /. z -> [Zeta], {{x, 0, 1}, {y, 0, 1}}]; But, unfortunately, I didn't get any solutions. Apr 6, 2019 at 3:40