Getting complex numbers results instead of expected real numbers

classic Classic list List threaded Threaded
5 messages Options
ant
Reply | Threaded
Open this post in threaded view
|

Getting complex numbers results instead of expected real numbers

ant
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!

OUTPUT:
************************************************************
octave:0>info_test
octave:1> info_vec =

    20.0000 +   0.0000i     9.6388 +   0.0000i
    30.0000 +   0.0000i    10.3678 -   0.0000i
    40.0000 +   0.0000i    11.2201 +   0.0000i
    50.0000 +   0.0000i    12.1703 +   0.0000i
    60.0000 +   0.0000i    13.2062 +   0.0000i
    70.0000 +   0.0000i    14.3192 -   0.0000i
    80.0000 +   0.0000i    15.5026 -   0.0000i
    90.0000 +   0.0000i    16.7513 -   0.0000i
   100.0000 +   0.0000i    18.0610 +   0.0000i
   110.0000 +   0.0000i    19.4281 +   0.0000i
   120.0000 +   0.0000i    20.8496 -   0.0000i
   130.0000 +   0.0000i    22.3228 -   0.0000i
   140.0000 +   0.0000i    23.8454 -   0.0000i
   150.0000 +   0.0000i    25.4154 +   0.0000i
   160.0000 +   0.0000i    27.0310 -   0.0000i
   170.0000 +   0.0000i    28.6904 +   0.0000i
   180.0000 +   0.0000i    30.3922 +   0.0000i
   190.0000 +   0.0000i    32.1350 -   0.0000i
   200.0000 +   0.0000i    33.9176 -   0.0000i
   210.0000 +   0.0000i    35.7389 -   0.0000i
   220.0000 +   0.0000i    37.5977 -   0.0000i

SCRIPT:
*****************************************************
global sigma r mu k0 alpha beta delta gama phistar c0
global l1p l1m l2p l2m

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));

gama = 10; phistar=0.6; alpha = 0.7; beta = 1/7.36936815669474; delta = 0.20;
c0_iter = linspace(20,220,21);

info_vec = [];
for c0 = c0_iter
    k0 = (beta*c0)^(1/alpha);
    function y1 = ClosureTrigger(xtau)
        global r mu sigma c0 gama alpha beta k0 phistar delta l1p l1m l2p l2m
        y1 = l2m*(r-mu)*(l1p-l1m)*(phistar*k0+(gama+delta*k0)/r-c0/r*(xtau/c0)^l2p)-((l2m-1)*(l1p-l1m)-l1m*(1+l1p-l2m)*(xtau/c0)^l1p)*xtau;
    endfunction
    function y2 = OptimalClosure()
        [y2] = fsolve(@ClosureTrigger,0);
    endfunction
    xstar = OptimalClosure();
    info_vec = [info_vec; c0 xstar];
endfor
info_vec
Reply | Threaded
Open this post in threaded view
|

Re: Getting complex numbers results instead of expected real numbers

Jordi Gutiérrez Hermoso-2
On 12 February 2012 09:19, ant <[hidden email]> 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!
[snip]
>        y1 =
> l2m*(r-mu)*(l1p-l1m)*(phistar*k0+(gama+delta*k0)/r-c0/r*(xtau/c0)^l2p)-((l2m-1)*(l1p-l1m)-l1m*(1+l1p-l2m)*(xtau/c0)^l1p)*xtau;

This expression can potentially raise negative numbers to powers. For
example, insert the following debugging lines in this function:

       if (iscomplex(y1))
         keyboard;
       endif

You will observe that this sometimes results in xtau/c0 being
negative, so you get a complex result with small imaginary part when
raising it to l1p.

If you want to avoid complex results when raising to powers, sometimes
the following functions can help: realpow, realsqrt, cbrt, nthroot.

HTH,
- Jordi G. H.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Getting complex numbers results instead of expected real numbers

CdeMills
In reply to this post by ant
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
ant
Reply | Threaded
Open this post in threaded view
|

Re: Getting complex numbers results instead of expected real numbers

ant
In reply to this post by ant
Thanks for your good suggestions, they are very useful for a beginner. However I'm still puzzled by the c0's showing as complex.

For example, since the c0 are all natural numbers (20, 30, 40, ..., 220) why I'm I getting complex numbers like 20 + 0000i,... in the end?

Letting c0 = 20 and changing the code to include one more line
        ...
        c0  #newline
    xstar = OptimalClosure()
    info_vec = [info_vec; c0 xstar];

OUTPUT:
c0 =  20
xstar = 9.6388e+000 + 6.0778e-054i
info_vec =
        20.0000 +  0.0000i    9.6388 +  0.0000i
       
My understanding is the following. Since info_vec has an element - xstar - that is complex, Octave "promotes" c0 to a complex number too, because now info_vec is a complex matrix. Is this correct?
Guess that, if so, this may cause quite a few surprises, and be difficult to debug if the complex numbers arise unexpectedly.
Reply | Threaded
Open this post in threaded view
|

Re: Getting complex numbers results instead of expected real numbers

CdeMills
Yes and no. In the construct:
info_vec = [info_vec; c0 xstar];

you concatenate a line to info_vec whose one element is complex, so the whole thing is promoted to complex.

Once again, the solutions is quite easy; in your function ClosureTrigger, just replace
y1 = l2m*(r-mu)*(l1p-l1m)*(phistar*k0+(gama+delta*k0)/r-c0/r*(xtau/c0)^l2p)-((l2m-1)*(l1p-l1m)-l1m*(1+l1p-l2m)*(xtau/c0)^l1p)*xtau;

by

y1 =real(l2m*(r-mu)*(l1p-l1m)*(phistar*k0+(gama+delta*k0)/r-c0/r*(xtau/c0)^l2p)-((l2m-1)*(l1p-l1m)-l1m*(1+l1p-l2m)*(xtau/c0)^l1p)*xtau);

This way, you drop the complex part which is just roundoff error. The whole call chain will then handle only real numbers; and the final results will be real too.

Regards

Pascal