Opened 10 years ago

Closed 2 years ago

#8423 closed enhancement (fixed)

Add code to plot Julia sets

Reported by: was Owned by: sdietzel
Priority: minor Milestone: sage-8.1
Component: dynamics Keywords: complexdynamics, gcoc2017
Cc: ncohen, brunellus, bhutz, atowsley Merged in:
Authors: Ben Barros Reviewers: Ben Hutz
Report Upstream: N/A Work issues:
Branch: 21e9be7 (Commits) Commit: 21e9be71bb46d349902c4f4aab0fe897c2c73b29
Dependencies: #23257 Stopgaps:

Description (last modified by bbarros)

Add code to sage for plotting Julia sets. This ticket is now a part of my Google Summer of Code project that is intended to improve the complex dynamics functionality in Sage. Please note that this ticket is now dependent upon #23257 (Plotting the Mandelbrot set in Sage). For more information on my project, please visit my blog:

Attachments (1)

trac-8423_part1.patch (2.6 KB) - added by was 10 years ago.
adds a "fractals" module to sage…

Download all attachments as: .zip

Change History (28)

Changed 10 years ago by was

adds a "fractals" module to sage...

comment:1 Changed 10 years ago by mhampton

Carl Witty and I wrote a function for integral curvature Apollonian gaskets: I have a Koch curve generator I can dig out as well.

comment:2 Changed 10 years ago by mhampton

Also, as a big "todo" list, this wikipedia page is pretty comprehensive:

comment:3 Changed 10 years ago by ncohen

  • Cc ncohen added

(cc: me)

comment:4 Changed 10 years ago by mhampton

Here's an iterated function system example (the fern):

def fern(x,y):
    An iterated function system whose orbit traces out
    a fern shape.
        x,y - numerical scalars
        a 2-component list of new x and y values.
    r = random()
    if r<.01:
        return [0,.16*y]
    elif r < .08:
        return [.2*x-.26*y, .23*x+.22*y+1.6]
    elif r < .15:
        return [-.15*x+.28*y, .26*x+.24*y+.44]
        return [.85*x+.04*y,-.04*x+.85*y+1.6]
def fern_orbit(n):
    Returns a trajectory of length n of the fern
    iterated function system.
        n - an integer, the length of the trajectory
        a list of 2-component lists
       sage: show(points(fern_orbit(10000),pointsize=1),axes=False)
    traj = [[0,0]]
    for i in range(n):
        nt = fern(traj[-1][0],traj[-1][1])
    return traj

comment:5 Changed 10 years ago by mhampton

The @interact example on the logistic bifurcation diagram might qualify as a fractal too:

comment:6 Changed 10 years ago by sdietzel

Actually, trac_8423_julia.patch is not correct. Don't use it.

comment:7 Changed 8 years ago by brunellus

  • Cc brunellus added

comment:8 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:9 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:10 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:11 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:12 Changed 3 years ago by kcrisman

Probably related: #23257

comment:13 Changed 2 years ago by bbarros

  • Authors set to Ben Barros
  • Cc bhutz atowsley added
  • Component changed from fractals to dynamics
  • Dependencies set to #23257
  • Description modified (diff)
  • Keywords complexdynamics gcoc2017 added
  • Milestone changed from sage-6.4 to sage-8.1
  • Status changed from new to needs_review
  • Summary changed from fractals: add code to plot julia sets to Add code to plot Julia sets

comment:14 Changed 2 years ago by bbarros

  • Branch set to u/bbarros/8423

comment:15 follow-up: Changed 2 years ago by bhutz

  • Commit set to 62fe4ece1ac4af6a011e6be9b257c8a23aa99696
  • Reviewers set to Ben Hutz
  • Status changed from needs_review to needs_work

I have only minor issues with this

  • I see that when you pass in a period you are giving a random root of the dynatomic. I'm not opposed to this, but it should definitely be documented. Especially since some of the roots are actually lower period points when the tail is not 0, for example julia_plot(period=[1,1]) returns a c value of a fixed point sometimes.
  • return_points - This seems to return the c value associated to a particular period. If I've understood this correctly, it is not needed.
  • In several places you say you are iterating c for the Julia set. You are not, you are iterating the complex point corresponding to that pixel for the fixed c value defining the map.

Last 10 new commits:

f74fac8Merge branch 'master' into 23257
e1a379123257: Changed to lazy_import in complex_dynamics/
1afa33223257: Fixed merge conflict and doctests
9b6738b23257: Cython, Python 3 improvements, Color widget added
a1d02e423257: Changed fast_mandel_plot to fast_mandelbrot_plot
8aae69c23257: Changed default for interact to False
f52fb3f23257: Added complex_dynamics to dynamics/index.rst
9ddbf6dMerge branch 'u/bbarros/23257' of git:// into 23257
848d10cMerge branch 'master' into 8423
62fe4ec8423: Added plotting Julia sets

comment:16 in reply to: ↑ 15 Changed 2 years ago by bbarros

The idea behind return_points was to provide a way for the user to find all of the Julia sets of a certain cycle structure instead of just getting a random one. For example, if I wanted to find all of Julia sets with cycle structure (2,3) I can run the following code:

c_values = julia_plot(period=[2,3], return_points=True)
for c in c_values:

I wasn't sure if this would be useful or not so if it isn't practical, I can remove it.

comment:17 Changed 2 years ago by bhutz

Yes, it would be nice to be able to generate those c values, but I don't think that belongs in the julia_plot function.

Having an example that loops through the roots of the dynatomic and plots each Julia set is probably a good example to add to the documentation so the user doesn't have to come up with their own way to do it.

comment:18 Changed 2 years ago by git

  • Commit changed from 62fe4ece1ac4af6a011e6be9b257c8a23aa99696 to 89e55c0913d157fa1a847beb41d4ad51ee60ed7b

Branch pushed to git repo; I updated commit sha1. New commits:

89e55c08423: Added julia_cycle_structure function

comment:19 Changed 2 years ago by bbarros

I separated the functionality of finding c values based on cycle structure into a separate function. Should this be something that users need to import before using? Also, to make sure that my documentation is correct, is the cycle structure being inputted for the point c under the map Q_c(z) = z * * 2 + c ?

comment:20 Changed 2 years ago by bhutz

Except that now you have the problem that this function will not always return correct answers (some the c values will have a different minimal period structure).

The right functionality here is something that takes a post-critical portrait and returns the maps with the specified family has that portrait (such as quadratic polynomials...). This is not a trivial construction in general and implementing this particular specialized function that sort of does it seems like not a good long term plan.

I'm still inclined to remove that function entirely and have an example that demonstrates how to do it as part of the julia_plot() documentation.

comment:21 Changed 2 years ago by bbarros

That sounds like a good idea. I'll remove the extra function and add the following example to the julia_plot documentation:

sage: period = [2,3] # not tested
....: R.<c> = QQ[]
....: P.<x,y> = ProjectiveSpace(R,1)
....: R = P.coordinate_ring()
....: H = End(P)
....: f = H([x^2+c*y^2,y^2])
....: L = f.dynatomic_polynomial(period).subs({x:0,y:1}).roots(ring=CC)
....: c_values = [k[0] for k in L]
....: for c in c_values:
....:     julia_plot(c)

Is that what you had in mind?

comment:22 Changed 2 years ago by git

  • Commit changed from 89e55c0913d157fa1a847beb41d4ad51ee60ed7b to deb85c845531a54627e9ae674c9776673b323e11

Branch pushed to git repo; I updated commit sha1. New commits:

deb85c88423: Removed julia_cycle_structure, added example

comment:23 Changed 2 years ago by bbarros

  • Status changed from needs_work to needs_review

comment:24 Changed 2 years ago by bhutz

Note that there is a merge conflict now.

Also, there are a few more minor things in the docs

  • misspelling - 'diplay'
  • You don't have all the instances of c and p right. It should be Q_c(p) and you are coloring p.
  • you need to wrap some of the kwds descriptions as those lines are too long
  • I'd also suggest saying that the period: with (formal) cycle structure to indicted that you are using the dynatomic polynomial
  • if period is not None I guess is ok, but I'm used to seeing if not period is None.

comment:25 Changed 2 years ago by git

  • Commit changed from deb85c845531a54627e9ae674c9776673b323e11 to 21e9be71bb46d349902c4f4aab0fe897c2c73b29

Branch pushed to git repo; I updated commit sha1. New commits:

789f496Merge branch 'master' into 8423
21e9be78423: Fixed documentation

comment:26 Changed 2 years ago by bhutz

  • Status changed from needs_review to positive_review

comment:27 Changed 2 years ago by vbraun

  • Branch changed from u/bbarros/8423 to 21e9be71bb46d349902c4f4aab0fe897c2c73b29
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.