Fractal to make and play with Python
Fractals are very good as teaching materials for math and programming. Introducing fractal figures that have been scratched while practicing.
I will not explain the theory.
Mandelbrot set
The entire code is here.
Speaking of fractals, it's the Mandelbrot set.
Definition $ z_0 = 0, z_ {n + 1} = z_n ^ 2 + A point on the complex plane where the complex sequence $ \ {z_n \} $ defined by C $ does not diverge at $ n \ to \ integer $ The set of C $ is called the Mandelbrot set.
Complex sequence here
Therefore, for each $ C $, $ \ {z_n \} $ is calculated for up to 100 items, and the convergence test is performed based on whether or not there is an absolute value exceeding 2. There are various coloring methods for the Mandelbrot set, but this time, we will color according to the number of iters required for the convergence test.
mandelbrot_set.py
solutions = np.array([2.7693, 0.11535-0.58974J, 0.11535+0.58974J])
def f(z, C):
return z*z + C
def convergence_judgment(C):
z = inital_z # 0.0
i = 0
while i < 100:
z = f(z, C)
if abs(z) > 2:
return i #Minimum of more than 2|z_n|Return n corresponding to.
i += 1
return i #If it does not diverge, return 100.
In the main function, write that the rectangle whose length is epsilon centered on the center coordinate center_position is subdivided into N × N grids, and the convergence test is performed in each grid.
mandelbrot_set.py
N = 1000
inital_z = 0
epsilon = 1 #<2
center_position = [0, 0]
def main():
grid = np.empty((N, N), dtype=np.uint8)
points_x = np.linspace(center_position[0]-epsilon, center_position[0]+epsilon, N)
points_y = np.linspace(center_position[1]-epsilon, center_position[1]+epsilon, N)
for n, ReC in tqdm(enumerate(points_x)):
for m, ImC in enumerate(points_y):
C = ReC + ImC*1J
grid[m, n] = convergence_judgment(C)
show(grid)
In this code, the calculation time was about 100 seconds with N = 1000. Click the image to enlarge.
 |
 |
 |
|---|---|---|
| (0,0) | (-0.53 -0.613) | (-0.76,0.067) |
It's out of the definition of the Mandelbrot set, but if you change the first term $ z_0 $, you can see the characteristic appearance again, which is fun! Also, it's fun to change the index 2 of z ^ 2 + C and see the pattern of z ^ n + C for any n!
Newton fractal
The entire code is here.
For example, a cubic equation for the complex number $ z $
Has three solutions. Since I am negligent, I can find a solution by entering the following in Mathics. 
Newton's method with any point on the complex plane as the starting point $ z_0 $
Consider updating $ z $ with to find the solution to equation (1). Of course, the left side of equation (1) is not monotonous with respect to $ z $, so we do not know if Newton's method will give us a solution. Newton's method results for a properly selected $ z_0 $
- Get to Solution 1
- Get to Solution 2
- Get to Solution 3
- Nowhere to reach (not converged by Newton's method)
You will get one of the following results. For each start point $ z_0 $ of Newton's method, the point $ z_0 $ on the complex plane is colored according to which result of 1 to 4 above is obtained.
newton_fractal.py
def f(z):
return z**3-3*z**2+z-1
def diff_f(z):
return 3*z**2-6*z+1
def find_nearest_solution_idx_and_residual(z):
nearest_solution_idx = abs(z-solutions).argmin()
residual = abs(z - solutions[nearest_solution_idx])
return nearest_solution_idx, residual
def solve_f(inital_z):
z = inital_z
for n in range(max_step):
z = z - f(z) / diff_f(z)
nearest_solution_idx, residual = find_nearest_solution_idx_and_residual(z)
if residual < residual_criteria: #convergence judgment
return nearest_solution_idx + 2*n/max_step
return len(solutions)
Then the following fractal will appear. The calculation time was 24 minutes with N = 2500.
 |
 |
 |
|---|
It's fun to play around with changing the equations!
[1] [Hiroyuki Inao Complex Dynamical Systems and Computability](https://www.math.kyoto-u.ac.jp/insei/?plugin=attach&pcmd=open&file=Inou.pdf&refer=KINOSAKI%20SEMINAR%202009% 2Fattach)
Sierpinski Gasket
The entire code is here
The gasket for the Sierpinski gasket is generated by the following algorithm.
- Determine the coordinates of the equilateral triangle (a1, a2, a3) appropriately.
- Connect the midpoints of each side of the equilateral triangle (a1, a2, a3).
- Of the four triangles that appeared in step 2, do 2 for the three triangles except the central triangle.
The source code below uses complex numbers for the coordinates for convenience, but this is not essential. A triangle is a three-element list, each element corresponding to the complex coordinates of the vertices of the triangle.
sierpinski_gasket_fractals.py
#Generate three triangles from one triangle.
def make_3triangles_from_triangle(triangle):
p0 = triangle[0]
p1 = triangle[1]
p2 = triangle[2]
new_p2 = (m*p0 + n*p1) / (m+n)
new_p0 = (m*p1 + n*p2) / (m+n)
new_p1 = (m*p2 + n*p0) / (m+n)
triangle0 = [p0, new_p1, new_p2]
triangle1 = [new_p1, p2, new_p0]
triangle2 = [new_p2, new_p0, p1]
return triangle0, triangle1, triangle2
def get_new_triangles(triangles):
new_triangles = []
for i, triangle in enumerate(triangles):
triangle0, triangle1, triangle2 = make_3triangles_from_triangle(triangles[i])
new_triangles.append(triangle0)
new_triangles.append(triangle1)
new_triangles.append(triangle2)
return new_triangles
Since it's a big deal, let's find the sum of the circumferences and the area. Theoretically, the circumference should converge to infinity and the area should converge to zero.
sierpinski_gasket_fractals.py
def get_circumference_and_area(triangles):
distance_between_two_point = abs(triangles[-1][0] - triangles[-1][1])
circumference = len(triangles) * distance_between_two_point
area = math.sqrt(3)/4.0 * distance_between_two_point ** 2
return circumference, area
When you run the program, you can see that the circumference increases and the area converges to 0.
depth:0 num_triangles:3 circumference:1.50000 area:0.10825 time:0.000010 [sec]
depth:1 num_triangles:9 circumference:2.25000 area:0.02706 time:0.000118 [sec]
depth:2 num_triangles:27 circumference:3.37500 area:0.00677 time:0.000203 [sec]
depth:3 num_triangles:81 circumference:5.06250 area:0.00169 time:0.000349 [sec]
depth:4 num_triangles:243 circumference:7.59375 area:0.00042 time:0.000826 [sec]
depth:5 num_triangles:729 circumference:11.39062 area:0.00011 time:0.001953 [sec]
depth:6 num_triangles:2187 circumference:17.08594 area:0.00003 time:0.006996 [sec]
depth:7 num_triangles:6561 circumference:25.62891 area:0.00001 time:0.044693 [sec]
depth:8 num_triangles:19683 circumference:38.44336 area:0.00000 time:0.117012 [sec]
depth:9 num_triangles:59049 circumference:57.66504 area:0.00000 time:0.377997 [sec]
drawing...
I made an animation because it was a big deal.
I changed the source code a little and made a strange Sierpinski gasket.
Again, the circumference diverges and the area converges to zero.
Dragon curve
The entire code is here.
As shown in the figure below, you can draw a curve picture called a dragon curve by adding images rotated 90 degrees around the end point of the line segment at each step.
If you rotate the point $ (x, y) $ around the point $ (C_x, C_y) $ by $ θ $, the point $ (x, y) $ will be
x_{new} = x\cos{\theta} - y\sin{\theta} + C_x - C_x\cos{\theta} + C_y\sin{\theta} \\
y_{new} = x\sin{\theta} + y\cos{\theta} + C_y - C_x\sin{\theta} - C_y\cos{\theta}
Go to. Make this easy
x_{new} = f_x(x,y,C_{x},C_{y}) \\
y_{new} = f_y(x,y,C_{x},C_{y})
I will write.
At a depth of $ j $, $ 2 ^ j $ new dots will appear. We will write the newly appearing $ 2 ^ j $ points as $ x_ {new} ^ {(k)} $, $ y_ {new} ^ {(k)} $. However, $ k = 1, 2,…, 2 ^ j $. Then, at the depth of $ j $, for each $ k $
x_{new}^{(k)} = f_x(x_{2^j-k},y_{2^j-k},x_{2^j},y_{2^j}) \\
y_{new}^{(k)} = f_y(x_{2^j-k},y_{2^j-k},x_{2^j},y_{2^j})
You can get the newly added coordinates by calculating.
DragonCurve.py
def get_rotate_coordinate(x, y, Cx, Cy, theta=Pi/2):
newx = x*math.cos(theta)-y*math.sin(theta)+Cx-Cx*math.cos(theta)+Cy*math.sin(theta)
newy = x*math.sin(theta)+y*math.cos(theta)+Cy-Cx*math.sin(theta)-Cy*math.cos(theta)
return newx, newy
def main():
# Initial Condition
xs = [1, 0]
ys = [1, 1]
for depth in range(max_depth):
for k in range(2**depth):
newx, newy = get_rotate_coordinate(xs[2**depth-k-1], ys[2**depth-k-1], Cx=xs[2**depth], Cy=ys[2**depth])
xs.append(newx)
ys.append(newy)
print("depth:{} length:{}".format(depth, len(xs)-1))
print("drawing..")
draw(xs, ys)
The figure obtained in this way is as follows.
 |
 |
 |
|---|---|---|
| depth 5 | depth 10 | depth 20 |
If you shift the rotation angle from $ \ theta = 90 ° $, the shape of the dragon curve will be continuously transformed, which is interesting!
Recommended Posts