ans =
solve() returning empty, does a solution exist?
Show older comments
Hey all. For a project I'm doing, I need to generate a random number from a custom probability density function. To do that, I am generating a random number using rand, and finding its corresponding y-value on the related cumulative distribution function. However, solve() and vpasolve() only seem to return answers for some of these random values. Looking at the graph of this function, however, it would seem that these values exist. The relevant code, as well as the cumulative distribution function Fx, are below. Is it possible that my cumulative distribution function is too complicated, or that a solution doesn't exist after all? Any help would be greatly appreciated! Also, I am new to Matlab, so apologies if my code contains rookie mistakes!
clear all;
%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 130;
stdev = 5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) = (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));
%Create PDF of x (fx)
fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);
%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
(x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));
%Example of solve that cannot be done
xval = solve(Fx == 0.0284, x);
% %Generate 10000 random x values via inverse CDF method, not working!
% randvals = rand(10000, 1);
% xvals = zeros(10000, 1, "sym");
% for i = 1:10000
% xvals(i) = solve(Fx == randvals(i), x);
% end
% yvals = fy(xvals);
% yvals = double(yvals);
2 Comments
For the second part of problem, use subs function to evaluate the function at different x values
clear all;
%Define the mean and standard deviation of the target normal distribution
%curve.
mu = 130;
stdev = 5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) = (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));
%Create PDF of x (fx)
fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);
%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
(x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));
%Example of solve that cannot be done
xval = solve(Fx == 0.0284,x);
% %Generate 10000 random x values via inverse CDF method, not working!
randvals = rand(10, 1);
for i = 1:10
xvals{i} = solve(Fx == randvals(i), x);
yvals{i} = double(subs(fy,x,xvals{i}));
end
yvals{:}
Before getting to the inverse sampling, we have to make sure the underlying pdf and cdf are correct.
Here is the original code, with slight modifications to enforce sym as much as possible
mu = sym(130);
stdev = sym(5)/sym(2);
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) = (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
normdist(p) = (1/(sqrt(2*sym(pi)*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));
%Create PDF of x (fx)
fx(x) = (normdist(fy) * abs(diff(fy)));
fx(x) cannot be a pdf because it doesn't integrate to unity
int(fx(x),x,-inf,inf)
vpa(ans)
Maybe fx(x) as defined should be divided by 3?
Given that fx(x) is defined as it is, there doesn't seem to be a reason for the piecewise manipulations. int works fine; it returns a piecewise expression on its own.
F1x(x) = int(fx(t),t,-inf,x);
vpa(F1x(x),5)
Here is original code
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);
%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
(x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));
vpa(Fx(x),5)
We see that Fx is not defined over two regions and returns NaN when evaluated in those regions
Fx([-20,20])
At least the limit of Fx(x) -> 1, as it should if it were a cdf
limit(Fx(x),x,inf)
Comparison shows quite a difference, even after normalizing F1
figure
fplot([Fx(x),F1x(x)/3],[-50,50])
Perhaps it would be helpful to explain in mathematical terms what the problem actually is.
Accepted Answer
John D'Errico
on 26 Jun 2026 at 22:48
Edited: John D'Errico
on 26 Jun 2026 at 22:53
Just because solve returns empty does not mean a solution does not exist. In fact, most relations you can write down have absolutely no analytical solution.
When solve does fail to find a nice pleasing analytical solution however, it then tries to pass the problem to vpasolve. And USUALLY, vpasolve can find a solution, if a solution exists for a nice, well-behaved real valued function. HOWEVER...
Again, however, vpasolve need not always succeed. It is just a numerical solver, and that means it uses a starting value. The default starting value as I recall, is a random number between 0 and 1. And sometimes I have seen a starting value in that interval will diverge. Let me see if I can come up with a simple example.
syms x real
f(x) = exp(x) + 5 - 70*exp(-(x-4)^2*750)
fplot(f,[-2,4.5])
grid on
That is a pretty simple function, one that has no analytical solution.
solve(f)
vpasolve(f)
As you can see, vpasolve also fails. I put a double solution in there, very near x==4. If I do start the solver in the vicinty of the known solution, it does succeed. The problem is, that solution is a difficult one to resolve, because it comes from a virtual spike in the function shape.
vpasolve(f,4)
In your case of course, your function is purely monotomic. And that means vpasolve SHOULD pretty much always suceed. So, looking more carefully at your problem, what happens?
mu = 130;
stdev = 5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) = (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));
%Create PDF of x (fx)
fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);
%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
(x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));
I do admit that Fx is virtually constant on the interval [0,1], which is again where the starting values for vpasolve will be. And that can potentially cause a problem, but not here, as long as you use vpasolve properly.
Fx(0) - 0.0284
Fx(1) - 0.0284
However, your problem is a subtly different one. Rather than trying to use vpasolve to find an exact solution to that equality statement, you want to SUBTRACT OFF 0.0284.
vpasolve(Fx(x) == 0.0284)
vpasolve(Fx(x) - 0.0284)
Indeed, that worked nicely. Why did it work, when the first call to vpasolve failed? Look at what MATLAB sees for these two expressions:
x0 = 6;
Fx(x0) == 0.0284
Fx(x0) - 0.0284
Do you see the equality statement will alsys return false for any number in the vicinity of the root? At the same time, the difference of Fx and a number returns a floating point number. Now vpasolve is a happy camper.
7 Comments
Aidan
on 27 Jun 2026 at 0:16
Walter Roberson
on 27 Jun 2026 at 0:33
vpasolve(Fx(x) - 0.0284 == 0) tries to solve for exact equality.
vpasolve(Fx(x) - 0.0284) looks for a numeric solution and is willing to accept a numeric approximation of zero.
Chances are that up to the maximum number of digits, that Fx(x) - 0.0284 produces an expression that is slightly below zero for candidate x, and is slightly greater than zero for the numerically adjacent x. vpasolve(Fx(x) - 0.0284) says that situation is good enough. vpasolve(Fx(x) - 0.0284 == 0) compares to exact 0 and says that the comparisons fail on both sides, so says there is no solution.
John D'Errico
on 27 Jun 2026 at 1:40
Edited: John D'Errico
on 27 Jun 2026 at 12:51
Again, look at the result returned by Fx(6) == 0.0284.
The result is either true or false. And so vpasolve will only ever be happy if the result is EXACTLY true.
The same problem arises when you write
Fx(x) - 0.0284 == 0
Again, this is either true or false. Remember this case?
x0 = 6;
Fx(6) == 0.0284
ans =
0.0084279305055892282445002708293545 == 0.0284
Look carefully at ans. What do you see? Not a number at all, but a statement of potential equality. Is the left hand side IDENTICALLY equal to the right hand side? Of course not. In fact, the two sides will be non-identical, and so the expression is false for essentially all but one real value of x. In this case, if we tried to evaluate that expression, it would return false, not a number, but a logical result, either true or false.
But vpasolve is a rootfinder. It wants to see a function that varies, and is smoothly differentiable. Everything is different for this:
Fx(x0) - 0.0284
ans =
-0.019972069494410771755499729170645
When you write the latter form, the result is a NUMBER. And when you do something like
vpasolve(fun)
then vpasolve looks for a zero of fun, whatever fun happens to be. This is not a question of a tolerance, but a question of what vpasolve sees, what you pass to vpasolve.
All of this is different when you use solve, since solve can look at an equality statement and understand it. That is, if we do this:
syms y
solve(y - 3 == 0)
ans =
3
solve can handle the equailty. But when solve is unable to find a zero, then it calls vpasolve. And vpasolve does not like the equality.
For this example, both versions work. Maybe the OP's definition of Fx is simply too difficult.
syms x
eqnLeft(x) = 200*sin(x);
eqnRight(x) = x^3 - 1;
s = vpasolve(eqnLeft(x)-eqnRight(x)==0.02, x)
s = vpasolve(eqnLeft(x)-eqnRight(x)-0.02, x)
After reading this thread multiple times, I'm still confused if there is a consensus as to whether or not these lines should or should not return (essentially, if not exactly) the same result
vpasolve(F(x) == c) % 1
vpasolve(F(x) - c == 0) % 2
vpasolve(F(x) - c) % 3
@Paul: VPASOLVE is a numeric solver. If the numerics were stored to infinite precision, then your three syntaxes would be equivalvent. But with finite precision, vpasolve(F(x)-c) is more likely to find a solution, for the reasons given.
This is fundamentally no different to any other numeric solver.
%vpasolve(F(x) == c) % 1
%vpasolve(F(x) - c == 0) % 2
%vpasolve(F(x) - c) % 3
The documentation states:
"If eqn is a symbolic expression (without the right side), the solver assumes that the right side is 0, and solves the equation eqn == 0."
That statement is unequivocal in that (2) and (3) should yield identical results. So I don't see why (3) is more likely to find a solution than (2).
But, this comment states that vpasolve treats (2) and (3) differently, which I don't see in my reading of the doc.
And given that vpasolve is a numeric solver, consider this comment that states vpasolve is looking for solutions that are "EXACTLY" or "IDENTICALLY" true for forms (1) and (2)? Surely, vpasolve has a tolerance on the accuracy of the solution (even if the doc doesn't clearly state what that tolelrance is).
Forms (1) and (2) are mathematicaly equivalent, differing by only a constant on either side of the equation, so if vpasolve finds a solution for one of them it should, in general, find a solution for the other. Depending on how vpasolve treats those two cases under the hood, I could kind of see why: a) the solutions returned from (1) and (2) might not be exactly the same, but they should be essentially the same, or b) maybe a corner case where one of those two forms fails due to the finitie precision. But I don't think we're talking about either of those situations in this thread.
Returning to the example of the OP:
mu = 130;
stdev = 5/2;
%Define the function whose input is being tested (fy), as well as the target
%normal distribution PDF (normdist).
syms x p t;
fy(x) = (-0.003867*x^3 + 0.01413*x^2 + 3.742*x + 101.9);
normdist(p) = (1/(sqrt(2*pi*(stdev)^2))*exp(-(p - mu)^2/(2*(stdev)^2)));
%Create PDF of x (fx)
fx(x) = (normdist(fy) * abs(diff(fy)));
%Temp fx function in terms of t for integrating
fxtemp(t) = fx(t);
%Create CDF of x (Fx)
%Determine integral bounds (if fy is not monotonic)
localextrema_xvals = vpasolve(diff(fy) == 0, x);
localextrema_yvals = fy(localextrema_xvals);
localmin = [localextrema_xvals(1) localextrema_yvals(1)];
localmax = [localextrema_xvals(2) localextrema_yvals(2)];
localmin_xints = vpasolve(fy == localmin(2), x, [-50 50]);
localmax_xints = vpasolve(fy == localmax(2), x, [-50, 50]);
interval1ub = localmax_xints(1);
interval3lb = localmin_xints;
%Set up function (piecewise if necessary)
interval1 = int(fxtemp, t, -Inf, interval1ub);
interval2 = int(fxtemp, t, localmin(1), localmax(1));
Fx(x) = piecewise(x < interval1ub, int(fxtemp, t, -Inf, x), ...
(x >= localmin(1)) & (x <= localmax(1)), interval1 + int(fxtemp, t, localmin(1), x),...
x > interval3lb, interval1 + interval2 + int(fxtemp, t, interval3lb, x));
Just to get an idea what Fx(x) looks like
vpa(Fx(x),5)
BTW, Fx(x) is incorrect as a CDF insofar as it is not defined for -34.784 <= x < -16.783 and for 19.219 < x <= 37.22. To wit:
Fx([-20,20])
figure
fplot(Fx(x),[0,10]) % region of interest
yline(0.0284)
Initial tries, vpasolve returns seemingly inconsistent results
xval1 = vpasolve(Fx == 0.0284,x)
xval2 = vpasolve(Fx - 0.0284 == 0,x)
xval3 = vpasolve(Fx - 0.0284,x)
From the plot, maybe x0 = 6 is a good initial guess that will return consistent results ....
xval1 = vpasolve(Fx == 0.0284,x,6)
xval2 = vpasolve(Fx - 0.0284 == 0,x,6)
xval3 = vpasolve(Fx - 0.0284,x,6)
but aparently not.
Specifying a range for the initial guess also doesn't help (does this also limit the search space? the doc needs to do better ...)
xval1 = vpasolve(Fx == 0.0284,x,[5,7])
xval2 = vpasolve(Fx - 0.0284 == 0,x,[5,7])
xval3 = vpasolve(Fx - 0.0284,x,[5,7])
From the plot, we know that the solution is somewhere around x = 6, so extract the portion of the piecewise Fx that applies
c = children((Fx(x)));
c22 = children(c{2,2});
c22lower = lhs(c22{1}); c22upper = rhs(c22{2});
Plot that section, but extend it to the right past the point where it's valid
figure;
fplot(c{2,1},double([c22lower,c22upper+10]))
Add lines to show where the solution should be found and the true upper bound on x that's applicable to c21
yline(0.0284);xline(double(c22upper))
Try all three forms with all three methods for applying the initial guess.
With no initial guess specified, vpasolve works fine and returns the exact same solution for all three forms, but finds the solution that's outside the region where c21 is valid.
xval1 = vpasolve(c{2,1} == 0.0284)
xval2 = vpasolve(c{2,1} - 0.0284 == 0)
xval3 = vpasolve(c{2,1} - 0.0284)
isAlways([xval1==xval2,xval1==xval3])
Specifying an initial point or range works fine and returns the exact same solution for all three forms
xval1 = vpasolve(c{2,1} == 0.0284,6)
xval2 = vpasolve(c{2,1} - 0.0284 == 0,6)
xval3 = vpasolve(c{2,1} - 0.0284,6)
isAlways([xval1==xval2,xval1==xval3])
xval1 = vpasolve(c{2,1} == 0.0284,[5,7])
xval2 = vpasolve(c{2,1} - 0.0284 == 0,[5,7])
xval3 = vpasolve(c{2,1} - 0.0284,[5,7])
isAlways([xval1==xval2,xval1==xval3])
These results suggest to me that vpasolve has some trouble due to the piecewise nature of Fx(x) (and maybe because Fx(x) is not defined over two intervals?), rather than anything to do with the specific form of the eqn input argument to vpasolve.
More Answers (0)
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!
