matlab fsolve求解非线性方程组为什么得到复数解

用matlab fsolve求解非线性方程组得到的是复数解，该方程组有实数解，而且发现求得的复数解的实部和答案基本吻合，这是为什么呢？下面是我建立的方程和fsolve求解的程序

function F = zw(X)
da= X(1);
dr= X(2);
ct= X(3);
dc=31;
Dw = 12.303;
A = 0.627; 
dm = 58;
a=35/180*pi;
Kn=3.0363*10^5;
ri=6.405;
Ri=dm/2+(ri-Dw/2)*cos(a);

    Fr=4500;Fa=3000;M=800000;
    p1 = 0:2*pi/13:24*pi/13;
    p2 = pi/13:2*pi/13:25*pi/13;
    F(1) = Fa-Kn*A^1.5*sum((sin(a) + da/A+Ri*(ct/A)*cos(p1)).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)+Kn*A^1.5*sum((sin(a) - da/A-Ri*(ct/A)*cos(p2)).*(((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*(ct/A)*cos(p2)).^2).^0.5);
    F(2) = Fr-Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p1)+0.5*dc*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)-Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p2)-0.5*dc*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) -da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5);
    F(3)=M-dm/2*Kn*A^1.5*sum((sin(a) + da/A+Ri*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)-dc/2*Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p1)+0.5*dc*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)+dm/2*Kn*A^1.5*sum((sin(a) - da/A-Ri*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*(ct/A)*cos(p2)).^2).^0.5)+dc/2*Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p2)-0.5*dc*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) -da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5);
  
end

h=optimset;
h.Display='off';
h.TolX=1e-6;
[jg,fval,exitflag]=fsolve('zw',[0,0,0],h) %求解

（1）问题分析
出现复数解应该是因为你的“zw”函数中有1.5次方和0.5次方，而fsolve搜索解的时候对应幂的底数可能是负的，这样对应的1.5次方和0.5次方就会出现复数，进而导致目标函数值为复数，最终导致得到的解为复数。
（2）问题的可能改进代码
考虑在zw中将最后的函数值取模，这样可以得到实数解

function F = zw_real(X)
da= X(1);
dr= X(2);
ct= X(3);
dc=31;
Dw = 12.303;
A = 0.627;
dm = 58;
a=35/180*pi;
Kn=3.0363*10^5;
ri=6.405;
Ri=dm/2+(ri-Dw/2)*cos(a);

Fr=4500;Fa=3000;M=800000;
p1 = 0:2*pi/13:24*pi/13;
p2 = pi/13:2*pi/13:25*pi/13;
F(1) = Fa-Kn*A^1.5*sum((sin(a) + da/A+Ri*(ct/A)*cos(p1)).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)+Kn*A^1.5*sum((sin(a) - da/A-Ri*(ct/A)*cos(p2)).*(((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*(ct/A)*cos(p2)).^2).^0.5);
F(2) = Fr-Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p1)+0.5*dc*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)-Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p2)-0.5*dc*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) -da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5);
F(3)=M-dm/2*Kn*A^1.5*sum((sin(a) + da/A+Ri*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)-dc/2*Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p1)+0.5*dc*(ct/A)*cos(p1)).*cos(p1).*(((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*ct*cos(p1)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5)+dm/2*Kn*A^1.5*sum((sin(a) - da/A-Ri*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) - da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*(ct/A)*cos(p2)).^2).^0.5)+dc/2*Kn*A^1.5*sum((cos(a) + (dr/A)*cos(p2)-0.5*dc*(ct/A)*cos(p2)).*cos(p2).*(((sin(a) -da/A-Ri*(ct/A)*cos(p2)).^2 + (cos(a) + dr* cos(p2)/A-0.5*dc*ct*cos(p2)/A).^2).^0.5 - 1).^1.5./((sin(a) + da/A+Ri*(ct/A)*cos(p1)).^2 + (cos(a) + dr* cos(p1)/A+0.5*dc*(ct/A)*cos(p1)).^2).^0.5);
F=abs(F);  %给函数返回值取模
end

clear;clc;
h=optimset;
h.Display='off';
h.TolX=1e-6;
format long;
[jg,fval,exitflag]=fsolve('zw_real',[0,0,0],h) %求解
format short;

（3）代码运行结果截图
这是将zw的返回值取模得到的结果，不知道能不能对上你的理论解，欢迎讨论

前提是你的方程组得有实数解，可以用assume函数或者syms后面加real来定义实数变量，然后用vpasolve函数求数值解，形如

syms x real
sol = vpasolve( [ ( x - 2 ) * ( x^2 + 1 ) == 0 ], [ x ] )

没实数解的话只能得到空集