How to get a random possible solution from 'solve' function when getting unknown parameters z1 and conditions
3 views (last 30 days)
Show older comments
Chen Yu Yuan
on 24 Feb 2024
Edited: John D'Errico
on 24 Feb 2024
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!!
0 Comments
Accepted Answer
John D'Errico
on 24 Feb 2024
Edited: John D'Errico
on 24 Feb 2024
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)
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])
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
That is, xrand will be some random value between 0 and 1.
vpasolve(subs(equations,x,xrand),[y,z])
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)
That is, xrand will be some random value between 0 and 1.
yzsol = vpasolve(subs(equations,x,xrand),[y,z])
yzsol.y
yzsol.z
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.
2 Comments
John D'Errico
on 24 Feb 2024
Edited: John D'Errico
on 24 Feb 2024
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))
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)
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.
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!