I m trying to solve non linear equations using fsolve function. Problem is I am getting the solution as imaginary. Is there any method to get only real valued solution?

function y=potFunc(pVect,eta,R,delC1,delC1p,delC2,delC2p)

tL=pVect(1);
tC=pVect(2);
d1=pVect(3);
d2=pVect(4);

s=sin(tL-tC);
f1=((d1*eta^2*sin(tL-tC)/(delC1p*cos(tL)))^2)^(1/3);
f2=((d2*eta^2*sin(tL-tC)/(delC2p*cos(tL)))^2)^(1/3);
y(1)=d1*s-delC1*cos(tL);
y(2)=d2*s-delC2*cos(tL);
y(3)=(d1*tan(tC)*(sin(tL)/sqrt(eta^2-sin(tL)^2)...
        -sqrt((f1-eta^2)/((eta^2-1)*f1)))-(R*sqrt((f1-eta^2)/((eta^2-1)))-tan(tC)))*1000;
y(4)=(d2*tan(tC)*(sin(tL)/sqrt(etya^2-sin(tL)^2)...
        -sqrt((f2-eta^2)/((eta^2-1)*f2)))-(R*sqrt((f2-eta^2)/((eta^2-1)))-tan(tC)))*1000;

%factor=norm(imag(pVect))^2; 

y=norm(y);
end

Output:

>> R=25;
>> eta=4/3;
>> delC1=0.98;
>> delC1p=0.87;
>> delC2=0.65;
>> delC2p=0.58;
>> d1=1.1;
>> d2=0.7;

>> pVect=[tL,tC,d1,d2];
>> f=@(pVect)potFunc(pVect,eta,R,delC1,delC1p,delC2,delC2p);
>> f(pVect)

ans =

      0.0968    0.1413   -0.0351   -0.1365

>> x=fsolve(f,pVect,optimset('MaxFunEvals',10000));

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 400 (the default value).

>> x

x =

     1.5717 - 0.0005i   1.5686 - 0.0750i   0.0060 + 0.0115i   0.0040 + 0.0077i

Answer 1

[Walter Roberson]

Use optimset to create an options structure that includes FunValCheck

  Check whether objective function values are valid. 'on' displays an error when the objective function returns a value that is complex, Inf, or NaN. The default, 'off', displays no error.

This will cause the solving to exit at the point where inputs that are purely real give a complex response from your potFunc(), and that will allow you to track down the bugs in your function.

fsolve() never spontaneously generates complex inputs, but once your function value goes complex on real inputs, the process of calculating the projected slopes and intercepts can end up with complex inputs because of the complex output from the function.