# Get random solution for equation with $n$ parameters and infinite solutions?

I'm trying to get a random solution for an equation with arbitrary number of parameters which has infinite solutions. Note that the answer must be in n dimensions where n is the number of parameters of the function.

My tries:

• The function can be something like x^2 + y^2 + z^2 == 1, in this case, the solution is a sphere and every point that satisfies that equality is a valid solution. This creates a region that consists of valid solutions, so if we generate a random point in that region we have a valid answer for this problem. That can be done using ybeltukov's answer's RegionDistribution:

implregion = ImplicitRegion[x^2 + y^2 + z^2 + t^2 == 1, {x, y, z, t}];
region = DiscretizeRegion[implregion];
pts = RandomVariate[RegionDistribution[region], 100]


but this raises an error saying DiscretizeRegion only takes arguments with less than 4 dimensions.

• I tried using NSolve to get all the points:

NSolve[x^2 + y^2 + z^2 == 1, {x,y,z}] (* I don't know what I was supposed to see here *)


but this gives a warning message that says:

 Warning: NSolve::infsolns: "Infinite solution set has dimension at least 1.
Returning intersection of solutions with (171802*x-113492*y-121484*z)/178835 == 1."


then it outputs some complex valued solutions, which are not valid (in this context):

{{x->0.252383 - 0.412144I, y->-1.17075 - 0.0795353I, z->-0.0214382 - 0.508549I},
{x->0.252383 + 0.412144I, y->-1.17075 + 0.0795353I, z->-0.0214382 + 0.508549I}}


As you can see, it gets those wanted infinite solutions, and it even selects a random solution, but it is always the same, and contains complex numbers when what I want is one random answer of that solution set which consists only of real values.

• I also tried with Reduce doing:

Reduce[x^2 + y^2 + z^2 == 1]


The result is:

(x == -Sqrt[1 - y^2 - z^2] || x == Sqrt[1 - y^2 - z^2])


and generating some points doing something like

Table[
{RandomChoice[{-1,1}] Sqrt[1 - y^2 - z^2], y, z}/.
{y-> RandomReal[{-1,1}],
z-> RandomReal[{-1,1}]},
{1000}]


gives me incorrect answers because x, y and z are not calculated "at the same time".

Edit:

• Using FindInstance distributes non uniformly the points, I used

FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, 1000]


and the result was • And the last method distributes non uniformly too, doing

Select[Table[{RandomChoice[{-1, 1}] Sqrt[1 - y^2 - z^2], y, z} /.
{y -> RandomReal[{-1, 1}],
z -> RandomReal[{-1, 1}]},
{1000}],
Im@#[] == 0 &]


gives me this is because the values are calculated depending on the $yz$ plane and in some places x changes faster than in another places.

• Have you considered FindInstance? Apr 15, 2017 at 18:58
• The results of your last approach, using Table, certainly are correct. What is your concern with them? Apr 15, 2017 at 19:00
• Would RandomPoint do what you're after? Apr 15, 2017 at 19:53
• @ChipHurst Yes! Apr 15, 2017 at 21:29
• Generating a uniform distribution even on the sphere is highly nontrivial, an ellipsoid is very hard, and you can forget about arbitrary hypersurfaces, so your question is much too general. Apr 15, 2017 at 23:43

FindInstance works quickly but seems to biased toward points that lie in the xy-plane.

SeedRandom;
With[{n = 10},
(FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, n] // N)[[All, All, 2]]]

{{0.790419, -0.108434, 0.602893}, {0.367265, 0.930116, 0.},
{-0.538922, -0.842356, 0.}, {-0.538922, 0.842356, 0.},
{0.55489, -0.831924, 0.}, {0.968064, 0.158416, -0.194311},
{-0.407186, -0.913345, 0.}, {0.978044, 0.0658436, 0.197724},
{0.946108, 0.323852, 0.}, {0.55489, 0.831924, 0.}}


So I recommend

SeedRandom;
With[{n = 10},
Module[{y, z},
z = RandomReal[{-1, 1}, n];
y = RandomReal[{-1 + #, 1 - #}] & /@ Abs[z];
MapThread[{Sqrt[1 - #1^2 - #2^2], #1, #2} &, {y, z}]]]

{{0.899547, -0.41092, -0.148189}, {0.972747, 0.0791143, -0.217954},
{0.903025, 0.301654, -0.305861}, {0.882229, 0.46164, -0.0925187},
{0.760152, 0.640032, 0.111927}, {0.701265, 0.574829, -0.421661},
{0.868464, 0.284055, -0.406304}, {0.714105, -0.38114, -0.587185},
{0.922918, 0.161119, -0.349661}, {0.321531, -0.0217459, 0.946649}}

• It is also biased, two peaks seem to have more points than the other two; check it out Apr 15, 2017 at 21:22

Region functionality in M11.2+ can handle Ellipsoid objects. For example:

pts = RandomPoint[
RegionBoundary @ Ellipsoid[{0,0,0,0},{1,2,3,4}],
10
]


{{0.702264, -0.506915, -1.00603, -2.29828}, {0.672747, 1.1246, 0.948138, -1.44966}, {-0.543024, 0.140708, -2.00356, -2.01652}, {0.652806, -1.37052, -0.93464, -0.339369}, {-0.38219, -0.873541, -2.1748, 1.48397}, {0.220089, -1.55027, 1.56574, -1.11951}, {-0.38265, 0.746885, 2.49541, -0.596283}, {0.52443, 1.36327, -1.48924, 0.471912}, {-0.212721, 1.43821, 1.50939, 1.71813}, {0.649815, -0.909925, 0.557857, -2.31921}}

Check that they all lie on the Ellipsoid:

Norm/@(pts.DiagonalMatrix[1/Range])


{1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}