8
$\begingroup$

I wish to plot all the points in 3 space that satisfy the following conditions.

$$ (x,y,z) = (a^2+b^2, ac+bd, c^2+d^2) $$ for $a,b,c,d$ in some range.

I have tried the following.

ParametricPlot3D[{a^2 + b^2, a*c + b*d, d^2 + c^2}, {a, -10, 
  10}, {b, -10, 10}, {c, -10, 10}, {d, -10, 10}]

However, it does not seem to work. Any suggestions would be much appreciated.

$\endgroup$
  • $\begingroup$ Doesn't that parameterization give you all points? $\endgroup$ – march Nov 22 '15 at 20:58
  • $\begingroup$ Are you sure there are spaces between parameters in your input, e.g. a c instead of ac? $\endgroup$ – Artes Nov 22 '15 at 21:01
  • 2
    $\begingroup$ You don't understand what you can do with ParametricPlot3D, you cannot visualize four parameter surfaces. $\endgroup$ – Artes Nov 22 '15 at 21:15
  • 1
    $\begingroup$ Related (nD -> 2D images): (13378), (24337), (26367) $\endgroup$ – Michael E2 Nov 23 '15 at 11:15
10
$\begingroup$

Method 1: unconstrained regions

You can easily do it with regions:

ℛ = ParametricRegion[{a^2 + b^2, a c + b d, d^2 + c^2}, {a, b, c, d}];

RegionPlot3D[ℛ, Axes -> True]

enter image description here

or

ineq = RegionMember[ℛ, {x, y, z}]
(* (x | y | z) ∈ 
  Reals && ((y == 0 && x >= 0 && z >= 0) || (z > 0 && -y^2 + x z >= 0)) *)

RegionPlot3D[ineq, {x, 0, 200}, {y, -200, 200}, {z, 0, 200}]

enter image description here

Method 2: numerical constrained regions

The original problem with restricted parameters is much more difficult (therefore, much more interesting). It is not possible to obtain exact region inequalities. However we can use several numerical tricks. Firstly, let us generate a huge set of points inside the region

region[a_, b_, c_, d_] := {a^2 + b^2, a c + b d, d^2 + c^2};
{{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} = {{-1, 1}, {-1, 1}, {-1, 1}, {-1, 1}};
{a, b, c, d} = {a1, b1, c1, 
    d1} + {a2 - a1, b2 - b1, c2 - c1, d2 - d1} Transpose@RandomReal[1, {2000000, 4}];
pts = Transpose@region[a, b, c, d];

For each point in space we can calculate the distance to the nearest point in the set

nf = Nearest@pts;
dist[x_?NumericQ, y_?NumericQ, z_?NumericQ] := 
  EuclideanDistance[First@nf@{x, y, z}, {x, y, z}];

We can plot a surface with small constant distance to the set (takes half a minute)

plot = RegionPlot3D[
  dist[x, y, z] < 0.1, {x, -0.1, 2.1}, {y, -2.1, 2.1}, {z, -0.1, 2.1}, 
  BoxRatios -> Automatic, PlotPoints -> {1, 2, 1} 20 + 1, Mesh -> {5, 10, 5}]

enter image description here

We got the approximate shape of the region. Fortunately, we can sufficiently improve the result. For each surface point we can calculate the nearest point in the region by FindArgMin. Initial values of parameters {a, b, c, d} are given by the corresponding parameters of the nearest point in the set pts. This procedure is a bit longer (takes several minutes)

nabcd = Nearest[pts -> Transpose@{a, b, c, d}];
nearest[x_, y_, z_] := 
  Block[{a, b, c, d}, 
     region @@ Quiet@
       FindArgMin[{Total[(region[a, b, c, d] - {x, y, z})^2], 
         a1 <= a <= a2 && b1 <= b <= b2 && c1 <= c <= c2 && 
          d1 <= d <= d2}, {{a, #1}, {b, #2}, {c, #3}, {d, #4}}, AccuracyGoal -> 10, 
        MaxIterations -> 20]] & @@ First@nabcd@{x, y, z};

plot /. GraphicsComplex[v_, 
   p_, ___] :> (GraphicsComplex[#, p, VertexNormals -> #2 - #] &[nearest @@@ v, v])

enter image description here

Now you see how thin is the region near edges (you even see mesh lines from the back surface).

Method 3: semi-analytical constrained regions

It is partially based on Michael's code. We can plot all 2-dimensional faces of the domain and analytically investigate the envelope. The envelope is given by $a d = b c$:

Resolve[
 Exists[{nx, ny, 
   nz}, {nx, ny, nz}.Grad[region[a, b, c, d], {a, b, c, d}] == {0, 0, 0, 
     0} && (nx != 0 || ny != 0 || nz != 0)], Reals]
(* (a == 0 && b == 0) || (a == 0 && c == 0) || (a != 0 && 
   c != 0 && -b c + a d == 0) || (c == 0 && d == 0) *)

It corresponds to two branches of $y = \pm \sqrt{xz} $. For each branch we can analytically calculate the region of possible $(x,z)$ pairs and plot it with Plot3D. The following function can do everything automatically (except some rare cases)

regionPlot[ranges_, meshstep_: 2] := 
  Module[{a, b, c, d, x, z, swapxy, f, vars, domain, domainfaces, xint, yint, zint, 
    envelope, faceimages, faceplots},
   f = {a^2 + b^2, a c + b d, d^2 + c^2};
   vars = {a, b, c, d};
   domain = Transpose@Prepend[Transpose@ranges, vars];
   domainfaces = Subsets[Range@Length@domain, {2}];
   swapxy = {{1, 0, 0}, {0, 0, 1}, {0, 1, 0}};
   envelope[s_] := With[{reg = ImplicitRegion[FullSimplify[
         Resolve[
          Exists[{a, b, c, d}, 
           x == f[[1]] && z == f[[3]] && a d == b c && s f[[2]] > 0 && 
            And @@ ((#2 < # < #3) & @@@ domain)], Reals], (x | z) ∈ 
          Reals], {x, z}]},
     If[Head@reg === EmptyRegion, {},
      GraphicsComplex[#.swapxy, ##2] /. (VertexNormals -> 
            v_) :> (VertexNormals -> -v.swapxy) & @@ 
       First@Quiet@
         Plot3D[s Sqrt[x z], {x, z} ∈ reg, 
          PlotRange -> All {{1, 1}, {1, 1}, {1, 1}}, Mesh -> None]]];

   faceimages = 
    Join @@ 
     Table[Prepend[Delete[domain, List /@ face], f /. Thread[vars[[face]] -> #]] & /@
        Tuples@ranges[[face]], {face, domainfaces}];

   faceplots = 
    First@ParametricPlot3D[##, PlotPoints -> 4, MaxRecursion -> 1, 
        Mesh -> (Range[##2, meshstep] & @@@ {##2}), BoundaryStyle -> Black] & @@@ 
     faceimages;

   Graphics3D[{faceplots, envelope /@ {1, -1}}, PlotRange -> All, 
    BoxRatios -> Automatic, Axes -> True]];

regionPlot[{{4, 10}, {4, 10}, {-10, 10}, {-10, 10}}]

enter image description here

Other possible results:

Do[Print@regionPlot[{2 Min@##, 2 + 2 Max@##} & @@@ 
    RandomInteger[{-5, 4}, {4, 2}]], {4}]

enter image description here enter image description here enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ The first way is neat, fast, but a bit rough. (+1) However, restricting the range on the parameters, as the OP seems to desire, seems harder; plotting this is taking forever: ParametricRegion[{a^2 + b^2, a c + b d, d^2 + c^2}, {{a, -10, 10}, {b, -10, 10}, {c, -10, 10}, {d, -10, 10}}]. $\endgroup$ – Michael E2 Nov 22 '15 at 23:25
  • $\begingroup$ Well, I've already upvoted :) I was just about to take a look at this again, but what you've got looks pretty good. $\endgroup$ – Michael E2 Nov 24 '15 at 1:18
8
$\begingroup$

Perhaps this is good enough:

ListPointPlot3D[
 Transpose[{a^2 + b^2, a c + b d, d^2 + c^2} /. 
   Thread[{a, b, c, d} -> Transpose@Tuples[Range[-10., 10., 1.], 4]]]]

Mathematica graphics


On the other hand, if we set $A=(a,b)$, $C=(c,d)$, $\alpha = ||A||$, $\gamma= ||C||$, then $$x = A\cdot A = \alpha^2,\ y = A \cdot C = \alpha\gamma\cos\theta, z = C \cdot C = \gamma^2\,,$$ where $\theta$ is the angle between $A$ and $C$. It follows that the point-set consists of all $(x,y,z)$ for which $y^2 \le x z$ with $x, z \ge 0$:

RegionPlot3D[y^2 <= x z, {x, 0, 200}, {y, -200, 200}, {z, 0, 200}]

Mathematica graphics


We can plot the envelope ($y^2 = xz$) and the images of the 2D faces of the domain:

f = {a^2 + b^2, a c + b d, d^2 + c^2};
domain = {{a, -10, 10}, {b, -10, 10}, {c, -10, 10}, {d, -10, 10}};
vars = domain[[All, 1]];
ranges = domain[[All, 2 ;; 3]];
domainfaces = Subsets[Range@Length@domain, {2}];

faceimages = DeleteDuplicatesBy[
    Sort[# /. {Thread[Variables[#] -> {s, t}], 
        Thread[Variables[#] -> {s, -t}], 
        Thread[Variables[#] -> {-s, t}], 
        Thread[Variables[#] -> {-s, -t}]}] &]@
   Flatten[Table[
     With[{sub = Thread[vars[[face]] -> Transpose@Tuples[ranges[[face]]]]},
      Thread[f /. sub]
      ],
     {face, domainfaces}],
    1];
faceplots = Function[face,
    ParametricPlot3D[face, 
     Evaluate[Sequence @@ Select[domain, MemberQ[Variables[face], First[#]] &]]]
    ] /@ faceimages;
envelope = ParametricPlot3D[#, {u, 0, 200}, {v, 0, 200}] & /@
    {{u, Sqrt[u v], v}, {u, -Sqrt[u v], v}};

Show[faceplots, envelope, PlotRange -> All]

Mathematica graphics

(Interestingly, the mesh-line bleed-through changes its appearance as the surface is rotated. Note that most unnecessary surfaces are deleted, but not all. Some internal surfaces are plotted.)

| improve this answer | |
$\endgroup$
  • $\begingroup$ Your last idea is very nice. Unfortunately, it is not general. For example, for {{a, 0, 10}, {b, 5, 10}, {c, -10, 10}, {d, -10, 10}} I get this (just to compare: img). $\endgroup$ – ybeltukov Nov 23 '15 at 17:12
  • $\begingroup$ @ybeltukov That's what I'd expect. You'd have to compute the part of the envelope that is included and treat the boundaries more generally than I did. (The assumption of equal ranges symmetric about zero is easily fixed, except that in full generality, it causes my kernel to crash frequently when done all in one command.) $\endgroup$ – Michael E2 Nov 23 '15 at 18:19
  • $\begingroup$ Michael E2, I generalized your approach a bit and put the rewritten code to my answer. I hope you don't mind, it can be interesting for you. $\endgroup$ – ybeltukov Nov 24 '15 at 1:13
  • $\begingroup$ @ybeltukov Yeah, I just saw it while you were writing this....Looks good. $\endgroup$ – Michael E2 Nov 24 '15 at 1:19
  • $\begingroup$ @Michael E2 Hi, How to change this code ListPointPlot3D[ Transpose[{a^2 + b^2, a c + b d, d^2 + c^2} /. Thread[{a, b, c, d} -> Transpose@Tuples[Range[-10., 10., 1.], 4]]]] to make a in range [1,2], b in range[1,3]... $\endgroup$ – Qi Zhong Nov 2 '16 at 2:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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