Fractal Related Discussion > Noob's Corner

(Question) Fractal in 3D associated with z^3-1

(1/2) > >>

r-lnr:
Hello guys,
Recently, I've been interested in doing 3D Fractals. After several attempts to code such a fractal, I come for help.
I have tried to use the following method : you have a fixe point P=(x,y,z) with z>0, for points of the form A=(x,y,w=0), you calculate w'=-k/u^3 for some k (e.g. 0.75), and then, if -k/u^3 is not near enough the altitude w you shift A on the line (AP) (and the direction, depending on w>w' or w'>w) and you iterat the process until w and w' are near enough.
To accelerate the computing of the points, I have used this reference http://www.fractalforums.com/mathematics/convergent-distance-estimation-t1566/ instead of using a constant shift distance.
However, I'm not familiar with fractals and thus I am wondering :

-How much iteration should the process of getting from A=(x,y,w=0) a point near enough the surface take ? I have been stopping the calculation at more than 1000 iterations with a precision of 0.0001 for w-w'. Do you think I should lower (at least in a first time) the precision ?

- For the distance estimation, it seems the relation is K=ln(|zn+1-zn|)|zn+1-zn|/|dzn+1-dzn| with dzn+1=dznP(zn)P''(zn)/(P'(zn)^2) and zn+1 and zn are near enough. However, i don't understand what value of dz0 should be taken (right now i have taken z1-z0), especially since it matters a lot in the estimation of the distance and so in my algorithm. Also should the shift be done of exactly the distance, or some constante times the distance (e.g. 0.9*distance) ?

-Do you know any cool ressources to learn that I could use that use this method ?

youhn:
Real fractals (as defined in maths) require an infinite number of iterations, just like integrals and derivatives require infinity for full precision. Fractals are the result of functions, and not all functions are analytical. Therefore some things are basically impossible to fully compute.

Luckily our time is finite, which makes things more pragmatic. Since "real" fractals do not exist in reality, we can only use approximations. Effectiveness is relative to purpose and desire. You don't need more exactness if the result is good enough.

Some resource for the 2D case:
https://en.wikipedia.org/wiki/Newton_fractal
https://mathematica.stackexchange.com/questions/149164/basins-of-attraction-when-using-newtons-method-to-solve-x3-1-0
http://paulbourke.net/fractals/newtonraphson/

Did you try any 2D fractals, before you started working on 3D fractals? Can you share any of your work, code or pretty pictures?

r-lnr:
I have done some 2D fractals, but it is a lot simplier to code, that's why I was interested in 3D fractals (and since what I did is not working, I was wondering whether there was any mistake or anything I did not understand in the process I tried to implement. Here is some fractal i obtained :

youhn:
Can you share the code from both the working 2D and the current 3D case for comparison?

r-lnr:
Here are both codes for 2D and 3D as asked

--- Code: ---from math import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LightSource

f= lambda z: z**3-1
fp= lambda z: 3*z**2
fpp=lambda z: 6*z
z1=1/2*1j*sqrt(3)-1/2
z2=-1/2*1j*sqrt(3)-1/2
z3=1
nx=50
ny=50
xmin=-3
xmax=3
ymin=-3
ymax=3
import numpy as np
x=np.linspace(xmin, xmax, nx)
y=np.linspace(ymin, ymax, ny)
x,y=np.meshgrid(x,y)
T=x+1j*y

nbitere=100
for n in range (nbitere) :
T=T-(f(T)/fp(T))

eps=0.01
T1=(abs(T-z1)<eps)
T2=(abs(T-z2)<eps)
T3=(abs(T-z3)<eps)

T=T1+2*T2+3*T3
import matplotlib.pyplot as plt
colo=np.array(['k','r','g','b'])
x= np.asarray(x).ravel()
y= np.asarray(y).ravel()
T= np.asarray(T).ravel()
couleur=colo[T]
plt.scatter(x,y, c=couleur)
plt.show()

##

PDV=0+1j*0
h=1
sup=1000
def temps(eps,z) :
t=1
a=z
d=-f(a)/(fp(a))
while t<sup and abs(a-z1)>eps and abs(a-z2)>eps and abs(a-z3)>eps :
t+=1
d=d*f(a)*fpp(a)/(fp(a)**2)
a=a-(f(a)/fp(a))
dis=abs(log(abs(f(a)/fp(a))))*abs(f(a)/fp(a))/(abs(d-d*f(a)*fpp(a)/(fp(a)**2)))
return(t,dis)

def doublew(u) :
return(-0.75/u**3)

xp=np.linspace(xmin, xmax, nx)
yp=np.linspace(ymin, ymax, ny)
xp,yp=np.meshgrid(xp,yp)
a=np.linspace(0,0,nx)
b=np.linspace(0,0,ny)
Z,Zp=np.meshgrid(a,b)
U=xp+1j*yp

eps2=0.01
C=0.9

for i in range(nx) :
j=0
print(i)
for j in range(ny) :
b=U[j][i]
w=0
t=0
c,dis=temps(eps,b)
while t<sup and abs(doublew(c)-w)>eps2 :
t+=1
if t==sup :
print ("pb")
if doublew(c)>w :
b=b+(C*dis*(PDV-b)/(abs(PDV-b)))
w=w+(C*dis*(h-w)/(abs(PDV-b)))
else :
b=b-(C*dis*(PDV-b)/(abs(PDV-b)))
w=w-(C*dis*(h-w)/(abs(PDV-b)))
c,dis=temps(eps,b)
xp[j][i]=b.real
yp[j][i]=b.imag
Z[j][i]=w

fig = plt.figure()

ax.plot_surface(xp, yp, Z)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

--- End code ---

Here I have taken nx and ny quite small, because it is already taking some time and not really working on small cases. The first part of the code is quite standard (it's Newton method), the second is the one I described in the first post (if it's not clear, do not hesitate to ask).