'Gnuplot not plotting negative values of x**(1/3)

This is the correct graph of x^(1/3):

Normal graph of x^(1/3)

However, when trying to plot x**(1./3) on gnuplot, this is what I get:

Graph of x^(1/3) using gnuplot

How can I fix this?

Additionally, it is not the first time that gnuplot doesn't plot values around zero in other functions (when it is supposed to, of course: I'm not talking about asymptotes etc.). What can I do?



Solution 1:[1]

The reason why your graph is not drawn closer to zero is because there is a default sampling of 100, i.e. your function will be plotted with 100 points within your xrange. If you increase this number, this will get you closer to zero. Check help samples.

About the missing negative values, this seems to be a bit special in gnuplot. I guess x**b where b is a floating point number is special. You can test the examples:

print sqrt(4)
print 4**0.5
print sqrt(-4)
print (-4)**0.5
print 8**(1./3)
print (-8)**(1./3)

which will give:

2.0
2.0
{0.0, 2.0}
{1.22460635382238e-16, 2.0}   # rounding error {0.0, 2.0}
2.0
{1.0, 1.73205080756888}       # 1 out of 3 valid solutions, but not the expected -2

Values in {Re, Im} are imaginary numbers. First the real part then the imaginary part. In order to get your plot nevertheless you can try the following:

Code:

### cube root
reset session

set samples 500
set grid xtics, ytics

cuberoot(x) = sgn(x)*abs(x)**(1./3)

plot cuberoot(x) w l
### end of code

Result:

enter image description here

Addition:

I'll try to explain, but I'm not a mathematician. For the N-th root there are N solutions. Apparently, gnuplot takes one of them. Apparently, one with a positive real part, and if there are several, the one with the smallest positive imaginary part. I guess it is called "principal root". Check also this.

So, this will explain why

print (-8)**(1./3)
print (-8)**(1./9)  # 9th root

will return

{1.0, 1.73205}          # and not -2
{1.18393, 0.4309183}    # and not -1.25992

Code:

### roots
reset session

set size ratio -1
set xlabel "real part"
set xrange [-3:3]
set ylabel "imaginary part"
set yrange [-3:3]
set grid xtics, ytics

set angle degrees
NRootRe(x,N,i) = -sgn(x)*abs(x)**(1./N)*cos(360.*i/N - 180*sgn(x))
NRootIm(x,N,i) = -sgn(x)*abs(x)**(1./N)*sin(360.*i/N - 180*sgn(x))

x = -8
do for [N=3:9:2] {
    do for [i=1:N] {
        print sprintf('x=%g, N=%g, i=%g: {% 9g, % 9g}',x,N,i,NRootRe(x,N,i), NRootIm(x,N,i))
    }
    print ""
}

plot for [N=3:9:2] [i=1:N:1] '+' u (0):(0):(NRootRe(x,N,i)):(NRootIm(x,N,i)) w vec ti sprintf("x=%g, N=%g",x,N)
### end of code

Result:

x=-8, N=3, i=1: {        1,  -1.73205}
x=-8, N=3, i=2: {        1,   1.73205}
x=-8, N=3, i=3: {       -2,  7.34764e-16}

x=-8, N=5, i=1: {-0.468382,  -1.44153}
x=-8, N=5, i=2: {  1.22624, -0.890916}
x=-8, N=5, i=3: {  1.22624,  0.890916}
x=-8, N=5, i=4: {-0.468382,   1.44153}
x=-8, N=5, i=5: { -1.51572,  5.56847e-16}

x=-8, N=7, i=1: {-0.839155,  -1.05227}
x=-8, N=7, i=2: { 0.299491,  -1.31216}
x=-8, N=7, i=3: {  1.21261, -0.583964}
x=-8, N=7, i=4: {  1.21261,  0.583964}
x=-8, N=7, i=5: { 0.299491,   1.31216}
x=-8, N=7, i=6: {-0.839155,   1.05227}
x=-8, N=7, i=7: {  -1.3459,  4.94459e-16}

x=-8, N=9, i=1: {-0.965156, -0.809862}
x=-8, N=9, i=2: {-0.218783,  -1.24078}
x=-8, N=9, i=3: { 0.629961,  -1.09112}
x=-8, N=9, i=4: {  1.18394, -0.430918}
x=-8, N=9, i=5: {  1.18394,  0.430918}
x=-8, N=9, i=6: { 0.629961,   1.09112}
x=-8, N=9, i=7: {-0.218783,   1.24078}
x=-8, N=9, i=8: {-0.965156,  0.809862}
x=-8, N=9, i=9: { -1.25992,  4.62872e-16}

enter image description here

Solution 2:[2]

The reason why Gnuplot doesn't plot or calculate x**(1/3.) for negative numbers (in real numbers world) is because it treats power function as logarithm functions: x^n = exp(ln(x^n)) = exp(nln(x)), where ln(x) is the natural logarithm. Notice that x^n = exp(n ln(x)) has its domain in x > 0.

Of course, we can compute ln(x) for x<0 using complex numbers (as we can see when compute (-8)^(1/3.))

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1
Solution 2