## My introduction to Compile

12

2

I am reading Fractals from the Newton-Raphson method from Peter Young's page.

I've tried:

newtC = Compile[{{n, _Integer}, {z, _Complex}},
Arg[FixedPoint[# - (#^n - 1)/(n #^(n - 1)) &, N[z], 50]]/(2 Pi)]


And:

Timing[DensityPlot[newtC[3, x + I y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 300]]


Gave me this image in 3.96397:

Now, I tried:

newt[n_, z_] :=
Arg[FixedPoint[# - (#^n - 1)/(n #^(n - 1)) &, N[z], 50]]/(2 Pi)


Then, even when I eliminated the PlotPoints->300, the following code will not work:

Timing[DensityPlot[newt[3, x + I y], {x, -2, 2}, {y, -2, 2}]]


I can't even abort the run: I have to quit the kernel. I'm really surprised at this huge difference. Am I missing something?

The docs have a similar example.

– J. M.'s ennui – 2015-11-21T16:51:53.843

9

Indeed, the newt function as you wrote it freezes Mathematica at least for a minute or so (I aborted it afterwards without waiting to see if it would complete).

Instead, you can prevent any attempts at symbolically evaluating the newt function by DensityPlot by restricting it to numerical arguments only:

Clear[newt]
newt[n_?NumericQ, z_?NumericQ] :=
Arg[FixedPoint[# - (#^n - 1)/(n #^(n - 1)) &, N[z], 50]]/(2 Pi)

Timing[DensityPlot[newt[3, x + I y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 300]]


Your compiled version of this function has a nice speed edge over the non-compiled one: on my system it takes only 3.5 s to generate the same density plot using your newtC.

Thanks. Just so everyone knows, this isn't mine, but it comes from a great page: http://young.physics.ucsc.edu/115/.

– David – 2015-11-21T06:04:57.277

4@David if you want to make the attribution more visible and permanent, then you should add it to your question. – Yves Klett – 2015-11-21T08:27:37.523

@YvesKlett Good point, I've done it. – David – 2015-11-21T15:29:00.027

@MarcoB I tried your approach and it only took 11.5376. You are saying restricting to numerical arguments only, yet the original function newt[n_, z_] := Arg[FixedPoint[# - (#^n - 1)/(n #^(n - 1)) &, N[z], 50]]/(2 Pi), if I change FixedPoint to FixedPointList, Clear[newt] newt[n_, z_] := Arg[FixedPointList[# - (#^n - 1)/(n #^(n - 1)) &, N[z], 50]]/(2 Pi), does seem to produce numeric arguments only: newt[3, -0.75 - 0.5 I]' produces{-0.406416, -0.321304, -0.339475, -0.333965, -0.333335, -0.333333,
-0.333333, -0.333333}. And there is the N[z] in the function. What's happening?
– David – 2015-11-21T15:44:22.567

3@David, as noted, DensityPlot[] will attempt to symbolically analyze your function to try to plot it better, which is not appropriate here. Putting in the _?NumericQ check thwarts that automatic analysis, and tells DensityPlot[]` to just concentrate on plotting it. – J. M.'s ennui – 2015-11-21T16:51:10.877