complex root of single variable nonlinear equation

I am solving the following single variable nonlinear equation.

syms x 

p=0;  % Require all roots for p=0 to p=25 with 0.1 interval

S=sqrt((p+sqrt(p^2+4*x^2))/2);

U=sqrt((-p+sqrt(p^2+4*x^2))/2);

Eq=S*U^5+S^5*U+2*S^3*U^3*cos(S)*cosh(U)+S^2*U^2*(S^2-U^2)*sin(S)*sinh(U);

ht = matlabFunction(Eq);

x0=2.4674;  % Initial guess for p=0

x = fzero(ht,x0);

For p=0 to 20.5, the roots are real but for p>20.5 the roots are complex (this is not a problem made by me, it is a text book problem of stability). As fzero gives only real solutions, there is no problem in solving till p=20.5, but for p=20.6 to 25 it gives some furthest (suddenly jumping to larger value) real root. My question is, how to get those complex roots. Can anyone please help me ?

Answer 1

[John D'Errico]

Easy enough.

syms x
p = 21;
S=sqrt((p+sqrt(p^2+4*x^2))/2);
U=sqrt((-p+sqrt(p^2+4*x^2))/2);
Eq=S*U^5+S^5*U+2*S^3*U^3*cos(S)*cosh(U)+S^2*U^2*(S^2-U^2)*sin(S)*sinh(U);

So solve, forcing vpasolve to look in the complex plane, and in a reasonable area.

r1 = vpasolve(Eq == 0,x,10+i)
r1 =
10.93735156115395089653090093632 + 2.1249406184976881586789820769267i

It appears the other root will be a complex conjugate of the first. But just to see if that is true, we can use a root killer:

vpasolve(Eq/(x-r1) == 0,x,10+i)
ans =
10.93735156115395089653090093632 - 2.1249406184976881586789820769267i

which finds the other root in that area. Did it work?

subs(Eq,x,r1)
ans =
0.00000000000000000000000000000099393633311179516308994541658166 + 0.00000000000000000000000000000062779071961032320818382550999783i

subs(Eq,x,r1')
ans =
0.00000000000000000000000000000099393633311179516308994541658166 - 0.00000000000000000000000000000062779071961032320818382550999783i

Yep.

[Kamal Bera]

It's working. Thanks.

Answer 2

[Matt J]

You would probably have to use fsolve instead

getx=@(z) complex(z(1),z(2)) ;

fun=@(z) [real(ht( getx(z)  )),    imag(ht(   getx(z)        ))];

    z0=[real(x0), imag(x0) ];
  z=fsolve(fun, z0);

x=getx(z),

[Matt J]

See also, Complex Numbers in Optimization Toolbox Solvers as well as this thread.