Finding real roots of polynomials

I want to find only real roots of the equation which is ;

4*sqrt((1-(z^2/f1^2))*(1-z^2))-(2-z^2)^2-(m*z^4*sqrt(1-z^2/f1^2)/sqrt(1-((z^2/f1^2)/y^2)))

I know that equation includes imaginary part and I do not want to see them. Moreover, my code fails and says that;

Error using fzero (line 242) Function values at interval endpoints must be finite and real.

Error in scholte (line 21) x=fzero(fun,x0)

Here is my code;

rho2=1000; %kg/m3
rho1=2700; %kg/m3
cl2=1481; %m/s
cl1=5919; %m/s
m=rho2/rho1;
y=cl2/cl1;
poi=0.25;

f1=(sqrt((1-2*poi)/(2*(1-poi))))^-1;

fun=@(z) 4*sqrt((1-(z^2/f1^2))*(1-z^2))-(2-z^2)^2-(m*z^4*sqrt(1-z^2/f1^2)/sqrt(1-((z^2/f1^2)/y^2)));

x0 = [1 10];
x=fzero(fun,x0)

[John D'Errico]

This is not a polynomial by the way. Not even a rational polynomial.

Answer 1

[Star Strider]

When in doubt, plot first:

fun=@(z) 4*sqrt((1-(z.^2/f1.^2)).*(1-z.^2))-(2-z.^2).^2-(m*z.^4.*sqrt(1-z.^2./f1.^2)./sqrt(1-((z.^2/f1.^2)./y.^2)));

z = linspace(1, 10);
fz = fun(z);
Refz = (imag(fz) == 0);                     % Indices Of Pure Real Values
fz_ext = [min(fz(Refz)) max(fz(Refz))];

figure(1)
plot(z(Refz), fz(Refz))
grid

There are no real roots on the interval [1 10].

[yusuf]

Thanks a lot. Then I should search new interval that includes real roots according to the plot, am I wrong ?

[Star Strider]

My pleasure.

I would plot over a new interval. When I run this:

[rts, fval] = fzero(fun, eps)

I get:

rts =

     -10.5367e-009

fval =

     444.0892e-018

So there is a root near zero. There may be more roots over a wider interval on ‘z’. I would plot it over the interval [-10 10] (or whatever interval interests you) to see what it does.

Answer 2

[John D'Errico]

There is a clear and exact root at z == 0.

And we have f1=sqrt(3) on the above computations. So with symbolic z, we get...

pretty(4*sqrt((1-(z^2/f1^2))*(1-z^2))-(2-z^2)^2-(m*z^4*sqrt(1-z^2/f1^2)/sqrt(1-((z^2/f1^2)/y^2))))
                                                             /      2 \
                                                       4     |     z  |
      /          /  2     \ \                      10 z  sqrt| 1 - -- |
      |   2      | z      | |     2     2                    \      3 /
4 sqrt| (z  - 1) | -- - 1 | | - (z  - 2)  - ----------------------------------
      \          \  3     / /                      /                       2 \
                                                   |     2251799813685248 z  |
                                            27 sqrt| 1 - ------------------- |
                                                   \       422926083573117   /

Kind of messy, but we can plot it to give some idea of what is happening.

ezplot(fun)
grid on

Hmm. That pretty much suggests if a root exists, then it is near zero. And vpasolve finds no other root besides zero.

vpasolve(fun)
ans =
0

That might suggest your efforts are best used to prove that no real root exists. Or is there one?

It seems best to first do a more detailed plot, in those regions where ezplot showed some strange behavior.

Ah, so there does appear to be a root, a bit under z=0.6. It turns out though, that this is just ezplot screwing up. In fact, above z=0.43 or so, it turns out that this function is always complex valued, but then it again shows real values above z=sqrt(3).

vpa(subs(fun,.4:.1:1),6)
ans =
[ 0.157388, 0.254125 + 0.0385171*i, 0.312266 + 0.0470278*i, 0.332791 + 0.0641264*i, 0.279062 + 0.0867165*i, 0.073598 + 0.114071*i, - 1.0 + 0.145422*i]

And finally, we can then find the second positive root.

ezplot(fun,[.42,.44])
grid on

vpasolve(fun,[.42 .44])
ans =
0.43256929327558654702065910622332

With a little effort, we can also show that for z just above sqrt(3), there is no real root to be found. My guess is, for larger values of z, it will be trivial to show the behavior is dominated by the higher powers of z, so that no other positive roots can possibly exist.

You could do similar investigations for negative z.

Answer 3

[yusuf]

Thanks a lot John. This would be solve my problem and give correct direction for this problem

[John D'Errico]

:)