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:

enter image description here

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?

David

Posted 2015-11-21T05:40:03.163

Reputation: 14 143

The docs have a similar example.

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

Answers

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]]

Mathematica graphics

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.

MarcoB

Posted 2015-11-21T05:40:03.163

Reputation: 53 573

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