Application Center - Maplesoft
  Submit your work
Application Center Applications Exponential map fractal viewer Preview

App Preview:

Exponential map fractal viewer

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



NULL

> 

NULL

> 

NULL

> 

NULL

> 

NULL

> 

NULL

> 

NULL

> 

restart

> 

with(plots):

> 

> 

Ox := 0;

> 

M := 1;

> 

w := 2/M:

> 

R := 400:

> 

lim := -1:

> 

Niter := 1000:

> 

Jiter := 300:

> 

> 

W := LambertW:

> 

eps := 0.8e-1:

> 

phi := proc (z) options operator, arrow; exp(z/exp(z)) end proc:

> 

delta := -0.2341e-1:

> 

k := 1;

> 

p := 3;

> 

t := (1+delta)*exp((2*k*Pi*I)/p)

> 

c := evalf(phi(t)):

> 

g := proc (z) options operator, arrow; c^z end proc:

> 

bail := 0.1e31:

> 

h := proc (k, z) options operator, arrow; -W(k, -log(z))/log(z) end proc:

> 

> 

HW := proc () local y, n, c, s, p, sol, i, aprx, dy, dist, r, newr, oldr, fun, dfun; if nargs < 2 then ERROR("At least two arguments required") end if; n := args[-1]; y := args[-2]; c := [args[1 .. -3]]; if y = 0 then 0 else dist := infinity; fun := exp(z); for i to nargs-2 do fun := exp(c[-i]*fun) end do; fun := z*fun-y; dfun := diff(fun, z); s := series(fun, z, n); p := convert(s, polynom); sol := {fsolve(p = 0, z, complex)}; for i to nops(sol) do aprx := evalf(subs(z = op(i, sol), fun)); dy := evalf(abs(aprx)); if dy <= dist then r := op(i, sol); dist := dy end if end do; oldr := r; newr := r-evalf(subs(z = r, fun)/subs(z = r, dfun)); for i to Niter while eps < abs((oldr-newr)/oldr) do oldr := newr; newr := newr-evalf(subs(z = newr, fun)/subs(z = newr, dfun)) end do; newr end if end proc:
> 

ppc := proc (c) local p, z0, at; global lim; p := -1; z0 := evalf(h(0, c)); at := evalf(abs(subs(z = z0, diff((g@@1)(z), z)))); if at <= 1 then p := 1; lim := z0 end if; if p = -1 then z0 := evalf(HW(-log(c), log(c), 15)/log(c)); at := evalf(abs(subs(z = z0, diff((g@@2)(z), z)))); if at <= 1 then p := 2; lim := z0 end if end if; if p = -1 then z0 := evalf(HW(-log(c), log(c), log(c), 15)/log(c)); at := evalf(abs(subs(z = z0, diff((g@@3)(z), z)))); if at <= 1 then p := 3; lim := z0 end if end if; if p = -1 then z0 := evalf(HW(-log(c), log(c), log(c), log(c), 15)/log(c)); at := evalf(abs(subs(z = z0, diff((g@@4)(z), z)))); if at <= 1 then p := 4; lim := z0 end if end if; if p = -1 then z0 := evalf(HW(-log(c), log(c), log(c), log(c), log(c), 25)/log(c)); at := evalf(abs(subs(z = z0, diff((g@@5)(z), z)))); if at <= 1 then p := 5; lim := z0 end if end if; p end proc:
> 

> 

lambda := evalf(log(c));

> 

ppc(c), lim;

> 

> 

RST := complexplot(phi(exp(I*theta)), theta = 0 .. 2*Pi, scaling = CONSTRAINED):

> 

cp := complexplot([c], style = POINT, symbol = DIAMOND):

> 

display(RST, cp);

> 

> 

J := proc (x, y) local z, m; z := evalhf(x+I*y); m := 1; to Jiter while eps < abs(z-lim) and abs(z) < bail do z := evalhf(exp(lambda*z)); m := m+1 end do; ln(m) end proc:
> 

> 

Jhf := proc (x, y) options operator, arrow; evalhf(J(x, y)) end proc:

> 

> 

t1 := time():

> 

p1 := plot3d(Jhf, Ox-w .. Ox+w, Oy-w .. Oy+w, grid = [R, R], shading = zhue, lightmodel = none, style = PATCHNOGRID, orientation = [-90, 0], axes = BOX, labels = [x, y, iter]):
> 

display(p1);

> 

t2 := time():

> 

print("Julia Set for c=", c, " in time:", t2-t1);

> 

> 

> 


> 

> 

> 

> 


> 


> 


> 


> 

> 

NULL

NULL