Why results of the 'solve' are wrong?

p_max=715.685424949238;
p_min=415.685424949238;
z_max=128.051406871193;
z_min=43.1985931288072;
D1_max=1128.05140687119;
D1_min=1043.19859312881;
angle=45;
y0=0:1:200;
p0=y0/sind(angle);  
x0=p0*cosd(angle); 
syms a b    
q=sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max==0;  %equation to solve  
z_solve=solve(q,b); 
z_value=subs(z_solve,a,p0); 
z0=eval(z_value);   
q1=sqrt((p0-p_min).^2+(z0-z_min).^2)-D1_min-sqrt((p0-p_max).^2+(z0-z_max).^2)+D1_max;   %check the reults

[Steven Lord]

You have not defined a variable named angle so when I run this code it tries to call the angle function with no input arguments and that throws an error. I recommend choosing a different name for your variable and editing your question to include the value of that variable.

[ruiheng wu]

sorry about that,'angle' is 45.

Answer 1

[Walter Roberson]

p_max=715.685424949238;

(and a bunch of other numbers).

What does that input mean? Does it mean that p_max is exactly ? Does it mean that p_max is some number not precisely defined that is between 715.6854249492375 (inclusive) and 715.6854249492385 (exclusive) ? Or is it intentional that p_max is exactly

fprintf('%.999g\n', p_max)

715.6854249492380404262803494930267333984375

which is how MATLAB will interpret it?

q=sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max==0;  %equation to solve  

You have sqrt() -- the Euclidean distance between (a,b) and (p_min, z_min) . most of the time, sqrt() is going to lead to irrational numbers. So your equation is going to try to exactly balance irrational numbers

z_solve=solve(q,b); 

Remember that the purpose of solve() is to find exact closed form algebraic solutions . So the solution it comes up with is not necesarily going to be remotely valid if p_max is exactly or if it is intended to be between 715.6854249492375 (inclusive) and 715.6854249492385 (exclusive). If you did not intend p_max to be exactly that number ending in "375" then you either formulated the problem incorrectly or else you made a mistake in using solve()

solve() is very much a "Garbage In, Garbage Out" solver. It is happy to give you the exact solution to the problem you actually asked it to solve, even though that might not be the problem you intended to solve. For any non-linear problem, do not assume that just because your inputs are "close" to the real inputs, that solve() is going to give you a solution that is "close" to the "real" solution: the exact solution can be a very fine balance between components of the expression.

[ruiheng wu]

Thank you for your explanation. I am giving a special case here. Next I will give the code applied to all cases and ask for your help to check if the problem is still caused by the same reason.

clear
angle=45;   %The results are correct when the angles are 30 degrees and 90 degrees.
f=1000;
y_max=550;
y_min=250;
R=(y_max-y_min)/2;
y_avg=(y_max+y_min)/2;
x_avg=y_avg*cosd(angle)/sind(angle);
y0=0:1:200;
d=sqrt(x_avg^2+y_avg^2);
p_max=d+R;
p_min=d-R;
z_max=p_max^2/(4*F);
z_min=p_min^2/(4*F);
D1_min=sqrt(p_min.^2+(F-z_min).^2);
D1_max=sqrt(p_max.^2+(F-z_max).^2);
p0=y0/sind(angle);  
x0=p0*cosd(angle); 
syms a b    
q=sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max==0;  %equation to solve  
z_solve=solve(q,b); 
z_value=subs(z_solve,a,p0); 
z0=eval(z_value);   
q1=sqrt((p0-p_min).^2+(z0-z_min).^2)-D1_min-sqrt((p0-p_max).^2+(z0-z_max).^2)+D1_max;   %check the reults

[Walter Roberson]

Q = @(v) sym(v);

syms F

angle = Q(45);   %The results are correct when the angles are 30 degrees and 90 degrees.
f = Q(1000);
y_max = Q(550);
y_min = Q(250);
R=(y_max-y_min)/2;
y_avg=(y_max+y_min)/2;
x_avg=y_avg*cosd(angle)/sind(angle);
y0 = Q(0:1:200);
d=sqrt(x_avg^2+y_avg^2);
p_max=d+R;
p_min=d-R;
z_max=p_max^2/(4*F);
z_min=p_min^2/(4*F);
D1_min=sqrt(p_min.^2+(F-z_min).^2);
D1_max=sqrt(p_max.^2+(F-z_max).^2);
p0=y0/sind(angle);  
x0=p0*cosd(angle); 
syms a b    
q=sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max==0;  %equation to solve  
z_solve = solve(q,b)

z_solve = 

z_value=subs(z_solve,a,p0); 
z0 = subs(z_value);   
q1=sqrt((p0-p_min).^2+(z0-z_min).^2)-D1_min-sqrt((p0-p_max).^2+(z0-z_max).^2)+D1_max;   %check the reults

simplify(q1(1:5))

ans = 

[ruiheng wu]

Thanks for your answer, by defining constants as symbolic vectors does avoid computational errors.

[Walter Roberson]

you probably intended a specific value for F, but you did not post it so I substituted symbolic F.

[ruiheng wu]

I made a mistake，'f' should have been 'F'.

[Walter Roberson]

Q = @(v) sym(v);

angle = Q(45);   %The results are correct when the angles are 30 degrees and 90 degrees.
F = Q(1000);
y_max = Q(550);
y_min = Q(250);
R=(y_max-y_min)/2;
y_avg=(y_max+y_min)/2;
x_avg=y_avg*cosd(angle)/sind(angle);
y0 = Q(0:1:200);
d=sqrt(x_avg^2+y_avg^2);
p_max=d+R;
p_min=d-R;
z_max=p_max^2/(4*F);
z_min=p_min^2/(4*F);
D1_min=sqrt(p_min.^2+(F-z_min).^2);
D1_max=sqrt(p_max.^2+(F-z_max).^2);
p0=y0/sind(angle);  
x0=p0*cosd(angle); 
syms a b    
q=sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max==0;  %equation to solve  
z_solve = solve(q,b)

z_solve = 

z_value=subs(z_solve,a,p0); 
z0 = subs(z_value);   
q1=sqrt((p0-p_min).^2+(z0-z_min).^2)-D1_min-sqrt((p0-p_max).^2+(z0-z_max).^2)+D1_max;   %check the reults

simplify(q1(1:5))

ans = 

Answer 2

[Paul]

Hi ruiheng,

I would stay out of the double world to first make sure everything is working as expected.

p_max=715.685424949238;
p_min=415.685424949238;
z_max=128.051406871193;
z_min=43.1985931288072;
D1_max=1128.05140687119;
D1_min=1043.19859312881;
angle=45;
y0=0:1:200;
p0=y0/sind(angle);  
x0=p0*cosd(angle); 
syms a b    
q = sqrt((a-p_min)^2+(b-z_min)^2)-D1_min-sqrt((a-p_max)^2+(b-z_max)^2)+D1_max == 0;  %equation to solve  

Add ReturnConditions

z_solve =  solve(q,b,'ReturnConditions',true);
vpa(z_solve.b,5)

ans = 

vpa(z_solve.conditions,5)

ans = 

%z_value = subs(z_solve,a,p0);
%z0 = eval(z_value)
%z0 = z_value(:,5);
%q1 = sqrt((p0-p_min).^2+(z0-z_min).^2)-D1_min-sqrt((p0-p_max).^2+(z0-z_max).^2)+D1_max   %check the reults
%simplify(q1)

Only checking a few solutions, timed-out when trying to check all

for ii = 8:12, %1:numel(p0)
    % verify that the value of a meets the condition for the solution
    cond(ii) = simplify(subs(z_solve.conditions,a,p0(ii)));
    % check the value of the solution for the value of a
    check(ii) = simplify(subs(lhs(q),[a b],{p0(ii) subs(z_solve.b,a,p0(ii))}));
end
cond(8:12)

ans = 

check(8:12)

ans = 

[ruiheng wu]

Thank you for checking the code.This should be the cause of the problem.