3

I am trying to plot the Newton's Basins of Attraction for the polynomial z^3-1 using python. I am using Newton-Raphson Iterative method to plot the graph.

Till now, I am able to plot this: Newton Basin My Plot

What I want to generate is something like this: Required Plot

Could someone please suggest how can I do it?


Update 1

After including the initial points (which I had mistakenly omitted): Improved Plot


Update 2

Here is the code for the Newton-Raphson loop. Is there any error in the code that causes the slow operation?

def newtonraphson(z):
    if z == 0:
        return False
    z1 = 0
    z2 = z
    tolerance = 1.0e-12

    while True:
        z1 = z2
        z2 = z1 - function(z1)/derivative(z1)
        if abs((z2 - z1).a) < tolerance:
            break

    return z2
Pratanu Mandal
  • 597
  • 8
  • 23
  • 1
    More detail is needed on the problem. Can you show your existing code? Do you think there is a specific problem that people might be able to help you with? Or is it that you don't understand the underlying concepts (which is probably beyond the scope of Stack Overflow)? – Stuart Nov 18 '15 at 18:01
  • How are you choosing your starting points in the complex plane? Your plot looks good, you just need to sample more points, I would say – Conor Nov 18 '15 at 18:13
  • 1
    @Conor I feel the same that I need to sample more points. The problem is, when I try to sample more points, the CPU usage reaches 100% and my computer starts lagging and sometimes hangs. How long should it take to generate such a plot on an average pc? – Pratanu Mandal Nov 18 '15 at 18:18
  • @Stuart I think I have understood the underlying concepts more or less and as you can see, my plot looks somewhat like the expected plot. I am passing all points by a step of 0.1 (..., -0.2, 0.1, 0.0, 0.1, 0.2, ...) for both real and imaginary parts to the Newton-Raphson iteration function and coloring the iteration path as per the final root. If I decrease the step to 0.05 (more sample points) my PC starts lagging and hangs as mentioned in my previous comment. – Pratanu Mandal Nov 18 '15 at 18:25
  • 1
    It would help if you at least added the code for the Newton loop that you use. At this resolution, several minutes for the plot is barely acceptable, any longer hints to a programming error. Assuming that the domain is [-3..3] in both directions, this are 60x60 points with 4-15 iterations of the Newton loop, or 4 times as many for step size 0.05, this should be really fast. – Lutz Lehmann Nov 18 '15 at 19:07
  • 1
    What complex numbers implementation do you use? In 'cmath' I do not find the `z.a` method. – Lutz Lehmann Nov 18 '15 at 19:41
  • @LutzL I have created my own complex number implementation. – Pratanu Mandal Nov 18 '15 at 19:43
  • 1
    Similar to this one https://stackoverflow.com/q/12486828/3088138? If `z.a` stands for the absolute value, why the extra `abs`? The Newton loop looks good, the quality change in the running time should have other reasons, perhaps in the data structure used to remember the computed points. – Lutz Lehmann Nov 18 '15 at 19:50
  • @LutzL No, the z.a stands for the real part of z 'abs' is for taking absolute value of the difference – Pratanu Mandal Nov 18 '15 at 20:06
  • As LutzL said, the code for Newton-Raphson looks right (although it's odd that you only check for the real component not changing to break the loop - you should check for the absolute value of the complex number not to change). – jcanizales Nov 18 '15 at 20:18
  • Because the problem is triggered when you increase the number of points x4, I would think it lies in how you're handling those points (outside the Netwon-Raphson solver). E.g. if you're leaking some big amount of memory per-point, augmenting the number of points might make your computer paginate. – jcanizales Nov 18 '15 at 20:20
  • Then it would be advisable to implement an absolute value method for your `Complex` class and use that. Complex numbers may also vary in the imaginary direction. – Lutz Lehmann Nov 18 '15 at 20:36
  • @LutzL OK. I shall implement an absolute value method for my complex class to account for the imaginary part. – Pratanu Mandal Nov 18 '15 at 20:41
  • @jcanizales Thank you so much! You have again pointed me in the right direction. Found a small bug (incorrect rounding off) which caused the problem. Solved it now! – Pratanu Mandal Nov 18 '15 at 20:49

1 Answers1

1

The plot you're getting looks "about right," except that it's lacking more points and them being evenly distributed.

How are you choosing which points to plot and which not to?

In general, you'll want to solve your equation your whole range of xs and ys as initial conditions. To select all those initial xs and ys, you want a small-enough step for your graph to look pretty (but not small enough that you're wasting time calculating points that end up being part of the same pixel).

Then you'll want to iterate the Newton-Raphson method not for a fixed number of steps, but rather until your solution is close enough to z^3-1 == 0 (equivalently, until z^3-1 is small enough). This way you're sure to get a color for each point.

EDIT to answer your comment:

when I try to sample more points, the CPU usage reaches 100% and my computer starts lagging and sometimes hangs

You need to give control back to the UI thread every once in a while. This will make your plot slower, but it will ensure that it is responsive. Ideally, you would perform all calculations in a background thread, and only communicate with your UI thread for it to plot each pixel.

jcanizales
  • 133
  • 1
  • 2
  • 11
  • If I do that (give back control to UI thread) or perhaps use sleep() (which I have already tried), the plot is taking a lot of time. I am passing all points by a step of 0.1 which completes in about 5-6 seconds. If I decrease the step to 0.05, it does not complete even within an hour! – Pratanu Mandal Nov 18 '15 at 18:30
  • Also, I am using matplotlib to plot the graph. Hence, the plot is not shown until the entire process is complete. So basically, everything happens in the background until finished. – Pratanu Mandal Nov 18 '15 at 18:32
  • 1
    <> you're not plotting all of them, though. In which cases are you not plotting a point? – jcanizales Nov 18 '15 at 18:35
  • If your calculations are happening in a background thread, your UI shouldn't hang. Even if the plot takes a long time to show up, you still should be able to move/minimize/maximize your window, close it, etc. – jcanizales Nov 18 '15 at 18:38
  • 1
    Your comment cleared everything. I was indeed not plotting all of them. I was mistakenly omitting the initial guesses. After incorporating them, the plot now does seem more like Newton's Basins of Attraction. Thank you! – Pratanu Mandal Nov 18 '15 at 18:58