Opened 10 years ago

Closed 2 years ago

# Add code to plot Julia sets

Reported by: Owned by: was sdietzel minor sage-8.1 dynamics complexdynamics, gcoc2017 ncohen, brunellus, bhutz, atowsley Ben Barros Ben Hutz N/A 21e9be7 (Commits) 21e9be71bb46d349902c4f4aab0fe897c2c73b29 #23257

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

### 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: http://neutraldrifts.blogspot.com/2009/01/integral-apollonian-packings-with-sage.html 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: http://en.wikipedia.org/wiki/List_of_fractals_by_Hausdorff_dimension

(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.

INPUT:
x,y - numerical scalars

OUTPUT:
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]
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 2-component 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 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: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
• Component changed from fractals to dynamics
• Dependencies set to #23257
• Description modified (diff)
• 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: ↓ 16 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:

 ​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 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:
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 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:

 ​89e55c0 `8423: 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:

 ​deb85c8 `8423: 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:

 ​789f496 `Merge branch 'master' into 8423` ​21e9be7 `8423: 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.