How to get a random possible solution from 'solve' function when getting unknown parameters z1 and conditions

Here are my codes:

syms x y z;
eq1 = x*x + y*y + z*z == 0.14;
eq2 = x*y + y*z + z*x == 0.11;

equations = [eq1, eq2];

assume(x>=0);
assumeAlso(y>=0);
assumeAlso(z>=0);
assumeAlso(x<=1);
assumeAlso(y<=1);
assumeAlso(z<=1);

S=solve(equations,[x,y,z],'ReturnConditions',true); 
disp(S)

S will print:

x: [4×1 sym]

y: [4×1 sym]

z: [4×1 sym]

parameters: z1

conditions: [4×1 sym]

and S.x is like:

(- 75*z1^2 - 30*z1 - 2)^(1/2)/10 - z1/2 - 3/10

- z1/2 - (- 75*z1^2 - 30*z1 - 2)^(1/2)/10 - 3/10

(- 75*z1^2 + 30*z1 - 2)^(1/2)/10 - z1/2 + 3/10

3/10 - (- 75*z1^2 + 30*z1 - 2)^(1/2)/10 - z1/2

and S.conditions is like:

5*z1 + 3 <= (- 75*z1^2 - 30*z1 - 2)^(1/2) & -13 <= 5*z1 + (- 75*z1^2 - 30*z1 - 2)^(1/2) & (- 75*z1^2 - 30*z1 - 2)^(1/2) <= 5*z1 + 13 & z1/2 + 3/10 <= -(- 75*z1^2 - 30*z1 - 2)^(1/2)/10 & 0 <= z1 & z1 <= 1

5*z1 + 3 <= (- 75*z1^2 - 30*z1 - 2)^(1/2) & -13 <= 5*z1 + (- 75*z1^2 - 30*z1 - 2)^(1/2) & (- 75*z1^2 - 30*z1 - 2)^(1/2) <= 5*z1 + 13 & z1/2 + 3/10 <= -(- 75*z1^2 - 30*z1 - 2)^(1/2)/10 & 0 <= z1 & z1 <= 1

5*z1 <= (- 75*z1^2 + 30*z1 - 2)^(1/2) + 3 & -7 <= 5*z1 + (- 75*z1^2 + 30*z1 - 2)^(1/2) & (- 75*z1^2 + 30*z1 - 2)^(1/2) <= 5*z1 + 7 & z1/2 <= 3/10 - (- 75*z1^2 + 30*z1 - 2)^(1/2)/10 & 0 <= z1 & z1 <= 1

5*z1 <= (- 75*z1^2 + 30*z1 - 2)^(1/2) + 3 & -7 <= 5*z1 + (- 75*z1^2 + 30*z1 - 2)^(1/2) & (- 75*z1^2 + 30*z1 - 2)^(1/2) <= 5*z1 + 7 & z1/2 <= 3/10 - (- 75*z1^2 + 30*z1 - 2)^(1/2)/10 & 0 <= z1 & z1 <= 1

How can i get any possible solution from parameter z1 and its conditions?

For instance, one of the possible solution would be x=0.1 y=0.2 z=0.3

But I have no idea how to get a random solution.

I'm sorry that my English is poor.

Thanks for your helping!!

Answer 1

[John D'Errico]

Solve cannot return a random possible solution. Wanting code to do what it is not designed to do will never be sufficient.

If you start from a specific start point though, vpasolve will generate a solution, at least some of the time.

So you can provide a random start point to vpasolve. This will work, at least if it is possible to work.

The problem is, you have only 2 equations, and 3 variables. So I'm actually surprised that solve was able to give you any solution at all.

syms x y z;
eq1 = x*x + y*y + z*z == 0.14;
eq2 = x*y + y*z + z*x == 0.11;
equations = [eq1, eq2];
assume(x>=0);
assumeAlso(y>=0);
assumeAlso(z>=0);
assumeAlso(x<=1);
assumeAlso(y<=1);
assumeAlso(z<=1);
S=solve(equations,[x,y,z],'ReturnConditions',true)

S = struct with fields:
             x: [4×1 sym]
             y: [4×1 sym]
             z: [4×1 sym]
    parameters: z1
    conditions: [4×1 sym]

It does mean that vpasolve will fail completely, as this is an under-determined problem. vpasolve generally does not handle that class of problem with more unknowns than equations.

vpasolve(equations,[x,y,z])

ans = struct with fields:
    x: [0×1 sym]
    y: [0×1 sym]
    z: [0×1 sym]

If you want a random solution to be generated, you can choose a specific randome value for one of the variables. For example...

xrand = rand

xrand = 0.9908

That is, xrand will be some random value between 0 and 1.

vpasolve(subs(equations,x,xrand),[y,z])

ans = struct with fields:
    y: [0×1 sym]
    z: [0×1 sym]

Of course this will fail most of the time. Why? Look at equation 1.

eq1 = x*x + y*y + z*z == 0.14;

So if x is greater than sqrt(0.14), NO real solution can ever exist. In fact, each of the variables x,y,z would follow the same constraint. Instead, try this:

xrand = rand*sqrt(0.14)

xrand = 0.2416

That is, xrand will be some random value between 0 and 1.

yzsol = vpasolve(subs(equations,x,xrand),[y,z])

yzsol = struct with fields:
    y: [2×1 sym]
    z: [2×1 sym]

yzsol.y

ans = 

yzsol.z

ans = 

And we now see a pair of solutions, generated conditionally subject to the value of x.

Can you do better than that? Perhaps, if we recognize that equation 1 is the equation of the surface of a sphere with three independent variables. It is centered at the origin, at the point (0,0,0).

How about equation 2? That is not so obvious, but it will generally be a hyperbolic surface. Plotting these equations in the first octant, we see:

fimplicit3(eq1,[0 1 0 1 0 1])
hold on
fimplicit3(eq2,[0 1 0 1 0 1])
hold off
grid on
box on
view(76,13)

The solution to the system of two equations is where those two surfaces intersect. Do you see the solution locus will be roughly a circle living in the 3-d space of your variables (x,y,z)? I say roughly a circle, because I cannot trivially prove it is circular. It appears to be circular, at least by my eye.

Knowing that fact, can we derive a general solution? Sigh. Possibly. I'm not at all sure it is worth the time I would need to invest though, on what appears to be possibly just an example problem.

[Chen Yu Yuan]

It really helps, thanks a lot!

This has been bothering me for a long time. Thank you for helping me resolve it.

Hope you have a wonderful day!

[John D'Errico]

I'm happy it helps. I might conjecture from the picuture that the solution locus lies in a plane. And that it is possibly circular. Given that as an idea, it might be possible to prove that assertion, knowing what to look for.

First, a quick numerical foray would be in order. I find it is often a good idea to verify my intuitions about a problem.

syms x y z;
eq1 = x*x + y*y + z*z == 0.14;
eq2 = x*y + y*z + z*x == 0.11;
equations = [eq1, eq2];
xyz = zeros(0,3);
while size(xyz,1) < 100
    xrand = rand*sqrt(0.14);
    yzsol = vpasolve(subs(equations,x,xrand),[y,z]);
    yz = double([yzsol.y,yzsol.z]);
    yz(imag(yz(:,1)) ~= 0,:) = [];
    if ~isempty(yz)
        nsol = size(yz,1);
        xyz = [xyz;[repmat(xrand,nsol,1),yz]];
    end
end

I've generated 100 solutions there. First, do they lie in a plane? This is easy enough to test.

format long g
S = svd(xyz - mean(xyz,1))

S = 3×1

1.13871106047224 0.836680188626194 2.41834573934206e-16

So they do indeed lie exactly in a plane, to within floating point trash. What is the equation of the plane that contains those solutions?

[U,S,V] = svd(xyz - mean(xyz,1));
V(:,3)

ans = 3×1

-0.577350269189626 -0.577350269189626 -0.577350269189626

Interesting. They lie in the plane with normal vector [1 1 1]. We might also notice that all solutions had the propery that x+y+z==0.6.

Do you see that is the solutions lie in a plane, then the intersectino of a plane with the sphere of equation 1 predicts that all solutions MUST lie on a perfect circle? Anyway, this is about where one should go back, and with that idea in hand, prove my claim with some mathematical rigor. I have developed a small case of rigor mortis, so I'll stop here.