1

Background

I read here that newton method fails on function x^(1/3) when it's inital step is 1. I am tring to test it in julia jupyter notebook.

  1. I want to print a plot of function x^(1/3)

  2. then I want to run code

f = x->x^(1/3)
D(f) = x->ForwardDiff.derivative(f, float(x))
x = find_zero((f, D(f)),1, Roots.Newton(),verbose=true)

Problem: How to print chart of function x^(1/3) in range eg.(-1,1)

I tried

f = x->x^(1/3)
plot(f,-1,1)

I got

enter image description here

I changed code to

f = x->(x+0im)^(1/3)
plot(f,-1,1)

I got

enter image description here

I want my plot to look like a plot of x^(1/3) in google enter image description here

However I can not print more than a half of it

enter image description here

Adam
  • 2,726
  • 1
  • 9
  • 22
grzegorzs
  • 486
  • 5
  • 11

3 Answers3

1

That Google plot makes no sense to me. For x > 0 it's ok, but for negative values of x the correct result is complex, and the Google plot appears to be showing the negative of the absolute value, which is strange.

Below you can see the output from Matlab, which is less fussy about types than Julia. As you can see it does not agree with your plot.

Plot of x^(1/3)

From the plot you can see that positive x values give a real-valued answer, while negative x give a complex-valued answer. The reason Julia errors for negative inputs, is that they are very concerned with type stability. Having the output type of a function depend on the input value would cause a type instability, which harms performance. This is less of a concern for Matlab or Python, etc.

If you want a plot similar the above in Julia, you can define your function like this:

f = x -> sign(x) * abs(complex(x)^(1/3))

Edit: Actually, a better and faster version is

f = x -> sign(x) * abs(x)^(1/3)

Yeah, it looks awkward, but that's because you want a really strange plot, which imho makes no sense for the function x^(1/3).

Community
  • 1
  • 1
DNF
  • 11,584
  • 1
  • 26
  • 40
1

That's because x^(1/3) does not always return a real (as in numbers) result or the real cube root of x. For negative numbers, the exponentiation function with some powers like (1/3 or 1.254 and I suppose all non-integers) will return a Complex. For type-stability requirements in Julia, this operation applied to a negative Real gives a DomainError. This behavior is also noted in Frequently Asked Questions section of Julia manual.

julia> (-1)^(1/3)
ERROR: DomainError with -1.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

julia> Complex(-1)^(1/3)
0.5 + 0.8660254037844386im

Note that The behavior of returning a complex number for exponentiation of negative values is not really different than, say, MATLAB's behavior

>>> (-1)^(1/3)

ans =

0.5000 + 0.8660i

What you want, however, is to plot the real cube root.

You can go with

plot(x -> x < 0 ? -(-x)^(1//3) : x^(1//3), -1, 1)

to enforce real cube root or use the built-in cbrt function for that instead.

plot(cbrt, -1, 1)

It also has an alias .

plot(∛, -1, 1)
hckr
  • 5,456
  • 1
  • 20
  • 31
1

F(x) is an odd function, you just use [0 1] as input variable.

The plot on [-1 0] is deducted as follow enter image description here

The code is below

import numpy as np

import matplotlib.pyplot as plt

# Function f 
f = lambda x: x**(1/3)


fig, ax = plt.subplots()

x1 = np.linspace(0, 1, num = 100)
x2 = np.linspace(-1, 0, num = 100)
ax.plot(x1, f(x1))
ax.plot(x2, -f(x1[::-1]))


ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')
plt.show()

Plot

enter image description here

Community
  • 1
  • 1
Adam
  • 2,726
  • 1
  • 9
  • 22