A complex number in the course of Müller's method

Consider the script below. I am trying to find the roots of by using Müller's method with initial estimations

N=input('Enter the degree of the main polynomial ');
C=input('Enter the constant ');
for k=1:N
a(k)=input('Enter a coefficient ');
end
a=[C,a];
Pol1=flip(a);
u=input('Enter the first estimate ');
v=input('Enter the second estimate ');
w=input('Enter the third estimate ');
z=input('Enter number of iterations ');
for n=1:z
delta0=(polyval(Pol1, v)-polyval(Pol1, u))/(v-u);
delta1=(polyval(Pol1, w)-polyval(Pol1, v))/(w-v);
h0=v-u;
h1=w-v;
a=(delta1-delta0)/(h1+h0);
b=a*h1+delta1;
c=polyval(Pol1,w);
g1=-2*c/(b+sqrt(b^2-4*a*c));   
g2=-2*c/(b-sqrt(b^2-4*a*c));
if b+sqrt(b^2-4*a*c)>b-sqrt(b^2-4*a*c)
d=g1;
end
if b+sqrt(b^2-4*a*c)<b-sqrt(b^2-4*a*c)
d=g2;
end
x(n)=w+d;  
u=v;
v=w;
w=x(n);
end

When I run it, MATLAB returns 'Undefined function or variable 'd'. d is calculated to be a complex number. MATLAB returns d when I run or separately, but does not do for the code. Why is that so?

Answer 1

[Walter Roberson]

if b+sqrt(b^2-4*a*c)>b-sqrt(b^2-4*a*c)
d=g1;
end
if b+sqrt(b^2-4*a*c)<b-sqrt(b^2-4*a*c)
d=g2;
end

That code will not assign d under a few circumstances:

• b+sqrt(b^2-4*a*c) == b-sqrt(b^2-4*a*c) -- which would occur if sqrt(b^2-4*a*c) == 0
• one of the values is nan
• sqrt(b^2-4*a*c) is imaginary but b is real; in that case, the > operator would compare only the real parts, and the real parts would be equal

Note:

b+sqrt(b^2-4*a*c)<b-sqrt(b^2-4*a*c) implies b + sqrt(b^2-4*a*c) - b < b - sqrt(b^2-4*a*c) - b implies sqrt(b^2-4*a*c) < - sqrt(b^2-4*a*c) implies 2*sqrt(b^2-4*a*c) < 0 implies sqrt(b^2-4*a*c) < 0 . However, algebraically, sqrt() is defined as "principle square root", which is positive when the value is positive, and is 1i * sqrt(-(b^2-4*a*c)) when b^2-4*a*c < 0 -- so sqrt(b^2-4*a*c) < 0 cannot happen unless you define an ordering of imaginary values with respect to 0. The MATLAB < operator ignores the imaginary component... but then you would be comparing 0 to 0, which would not be < .

Which is to say that your second test can never succeed, and your first test is always true unless b^2-4*a*c == 0 (a ligitimate test) or b^2-4*a*c < 0 (because the imaginary component is ignored)

[Ali Kiral]

A thorough explanation. Thank you, Walter