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:  sage8.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 )
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: https://benbarros.wordpress.com
Attachments (1)
Change History (28)
Changed 10 years ago by
comment:1 Changed 10 years ago by
Carl Witty and I wrote a function for integral curvature Apollonian gaskets: http://neutraldrifts.blogspot.com/2009/01/integralapollonianpackingswithsage.html I have a Koch curve generator I can dig out as well.
comment:2 Changed 10 years ago by
Also, as a big "todo" list, this wikipedia page is pretty comprehensive: http://en.wikipedia.org/wiki/List_of_fractals_by_Hausdorff_dimension
comment:4 Changed 10 years ago by
Here's an iterated function system example (the fern):
def fern(x,y): """ An iterated function system whose orbit traces out a fern shape. INPUT: x,y  numerical scalars OUTPUT: a 2component 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] else: 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. INPUT: n  an integer, the length of the trajectory OUTPUT: a list of 2component lists EXAMPLE: 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]) traj.append(nt) return traj
comment:5 Changed 10 years ago by
The @interact example on the logistic bifurcation diagram might qualify as a fractal too:
http://wiki.sagemath.org/interact/dynsys#CythonizedLogisticOrbitMap
comment:6 Changed 10 years ago by
Actually, trac_8423_julia.patch is not correct. Don't use it.
comment:7 Changed 8 years ago by
 Cc brunellus added
comment:8 Changed 6 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:9 Changed 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:10 Changed 6 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:11 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:12 Changed 3 years ago by
Probably related: #23257
comment:13 Changed 2 years ago by
 Cc bhutz atowsley added
 Component changed from fractals to dynamics
 Dependencies set to #23257
 Description modified (diff)
 Keywords complexdynamics gcoc2017 added
 Milestone changed from sage6.4 to sage8.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
 Branch set to u/bbarros/8423
comment:15 followup: ↓ 16 Changed 2 years ago by
 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:
f74fac8  Merge branch 'master' into 23257

e1a3791  23257: Changed to lazy_import in complex_dynamics/all.py

1afa332  23257: Fixed merge conflict and doctests

9b6738b  23257: Cython, Python 3 improvements, Color widget added

a1d02e4  23257: Changed fast_mandel_plot to fast_mandelbrot_plot

8aae69c  23257: Changed default for interact to False

f52fb3f  23257: Added complex_dynamics to dynamics/index.rst

9ddbf6d  Merge branch 'u/bbarros/23257' of git://trac.sagemath.org/sage into 23257

848d10c  Merge branch 'master' into 8423

62fe4ec  8423: Added plotting Julia sets

comment:16 in reply to: ↑ 15 Changed 2 years ago by
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: julia_plot(c)
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
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
 Commit changed from 62fe4ece1ac4af6a011e6be9b257c8a23aa99696 to 89e55c0913d157fa1a847beb41d4ad51ee60ed7b
Branch pushed to git repo; I updated commit sha1. New commits:
89e55c0  8423: Added julia_cycle_structure function

comment:19 Changed 2 years ago by
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
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 postcritical 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
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
 Commit changed from 89e55c0913d157fa1a847beb41d4ad51ee60ed7b to deb85c845531a54627e9ae674c9776673b323e11
Branch pushed to git repo; I updated commit sha1. New commits:
deb85c8  8423: Removed julia_cycle_structure, added example

comment:23 Changed 2 years ago by
 Status changed from needs_work to needs_review
comment:24 Changed 2 years ago by
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
andp
right. It should beQ_c(p)
and you are coloringp
.
 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 seeingif not period is None
.
comment:25 Changed 2 years ago by
 Commit changed from deb85c845531a54627e9ae674c9776673b323e11 to 21e9be71bb46d349902c4f4aab0fe897c2c73b29
comment:26 Changed 2 years ago by
 Status changed from needs_review to positive_review
comment:27 Changed 2 years ago by
 Branch changed from u/bbarros/8423 to 21e9be71bb46d349902c4f4aab0fe897c2c73b29
 Resolution set to fixed
 Status changed from positive_review to closed
adds a "fractals" module to sage...