Re: Getting complex numbers results instead of expected real numbers

[CdeMills]

ant wrote
> 
> When running a script in Octave 3.4.3, I'm getting complex numbers of the
> form (a+0.0000i) or (a-0.0000i), instead of the expected real number a.
> Any help will be much appreciated!
> 
> sigma = 0.10; mu = -0.02; r = 0.10;
> l1p = -mu/sigma^2 -1/2 + sqrt((mu/sigma^2+1/2)^2+(2*(r-mu)/sigma^2));
> l1m = -mu/sigma^2 -1/2 - sqrt((mu/sigma^2+1/2)^2+(2*(r-mu)/sigma^2));
> l2p = -mu/sigma^2 +1/2 + sqrt((mu/sigma^2-1/2)^2+(2*r/sigma^2));
> l2m = -mu/sigma^2 +1/2 - sqrt((mu/sigma^2-1/2)^2+(2*r/sigma^2));
> 
> 

What I see is that sometimes the argument of sqrt should be zero, and due to
rounding error you get something like -eps. This generates the imaginary
part of value j*sqrt(eps).

What you can do:
1) there are a lot of repeat in your code. Do something like
tmp = (mu/sigma^2+1/2)^2+(2*(r-mu)/sigma^2));
idx_OK = abs(tmp) > eps; tmp(idx_OK) = sqrt(tmp); tmp(~idx_OK) = 0;
l1p = mu/sigma^2 -1/2 + tmp;
l1m = mu/sigma^2 -1/2 - tmp;

2) just throw away the imaginary part
l1p = real(-mu/sigma^2 -1/2 + sqrt((mu/sigma^2+1/2)^2+(2*(r-mu)/sigma^2)));

Regards

Pascal


--
View this message in context: 
http://octave.1599824.n4.nabble.com/Getting-complex-numbers-results-instead-of-expected-real-numbers-tp4381221p4381333.html
Sent from the Octave - General mailing list archive at Nabble.com.