0
$\begingroup$

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]
Improve this question
$\endgroup$
2
  • $\begingroup$ show you initial code of before solution! to make our life easier ! $\endgroup$
    – nufaie
    Apr 4, 2019 at 23:35
  • $\begingroup$ @Alrubaie sure thing! $\endgroup$
    – Lucas Cruz
    Apr 5, 2019 at 2:24

1 Answer 1

Reset to default
0
$\begingroup$

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]

enter image description here

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.
enter image description here

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]
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – Lucas Cruz
    Apr 5, 2019 at 13:30
  • $\begingroup$ 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. $\endgroup$
    – LouisB
    Apr 5, 2019 at 18:47
  • $\begingroup$ 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. $\endgroup$
    – Lucas Cruz
    Apr 6, 2019 at 3:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.