1
1
from typing import Tuple
2
+
2
3
import numpy as np
3
4
from numpy .polynomial import Polynomial
4
5
5
- def mandelbrot (width : int , height : int , * ,
6
- x : float = - 0.5 , y : float = 0 , zoom : int = 1 , max_iterations : int = 100 ) -> np .array :
6
+
7
+ def mandelbrot (
8
+ width : int ,
9
+ height : int ,
10
+ * ,
11
+ x : float = - 0.5 ,
12
+ y : float = 0 ,
13
+ zoom : int = 1 ,
14
+ max_iterations : int = 100
15
+ ) -> np .array :
7
16
"""
8
17
From https://www.learnpythonwithrune.org/numpy-compute-mandelbrot-set-by-vectorization/.
9
18
"""
10
19
# To make navigation easier we calculate these values
11
- x_width , y_height = 1.5 , 1.5 * height / width
12
- x_from , x_to = x - x_width / zoom , x + x_width / zoom
13
- y_from , y_to = y - y_height / zoom , y + y_height / zoom
20
+ x_width , y_height = 1.5 , 1.5 * height / width
21
+ x_from , x_to = x - x_width / zoom , x + x_width / zoom
22
+ y_from , y_to = y - y_height / zoom , y + y_height / zoom
14
23
15
24
# Here the actual algorithm starts
16
25
x = np .linspace (x_from , x_to , width ).reshape ((1 , width ))
17
26
y = np .linspace (y_from , y_to , height ).reshape ((height , 1 ))
18
- c = x + 1j * y
27
+ c = x + 1j * y
19
28
20
29
# Initialize z to all zero
21
30
z = np .zeros (c .shape , dtype = np .complex128 )
@@ -26,27 +35,38 @@ def mandelbrot(width: int, height: int, *,
26
35
# To keep track on which points did not converge so far
27
36
m = np .full (c .shape , True , dtype = bool )
28
37
for i in range (max_iterations ):
29
- z [m ] = z [m ]** 2 + c [m ]
30
- diverged = np .greater (np .abs (z ), 2 , out = np .full (c .shape , False ), where = m ) # Find diverging
31
- div_time [diverged ] = i # set the value of the diverged iteration number
32
- m [np .abs (z ) > 2 ] = False # to remember which have diverged
38
+ z [m ] = z [m ] ** 2 + c [m ]
39
+ diverged = np .greater (
40
+ np .abs (z ), 2 , out = np .full (c .shape , False ), where = m
41
+ ) # Find diverging
42
+ div_time [diverged ] = i # set the value of the diverged iteration number
43
+ m [np .abs (z ) > 2 ] = False # to remember which have diverged
33
44
34
45
return div_time
35
46
36
- def julia (width : int , height : int , * ,
37
- c : complex = - 0.4 + 0.6j , x : float = 0 , y : float = 0 , zoom : int = 1 , max_iterations : int = 100 ) -> np .array :
47
+
48
+ def julia (
49
+ width : int ,
50
+ height : int ,
51
+ * ,
52
+ c : complex = - 0.4 + 0.6j ,
53
+ x : float = 0 ,
54
+ y : float = 0 ,
55
+ zoom : int = 1 ,
56
+ max_iterations : int = 100
57
+ ) -> np .array :
38
58
"""
39
59
From https://www.learnpythonwithrune.org/numpy-calculate-the-julia-set-with-vectorization/.
40
60
"""
41
61
# To make navigation easier we calculate these values
42
- x_width , y_height = 1.5 , 1.5 * height / width
43
- x_from , x_to = x - x_width / zoom , x + x_width / zoom
44
- y_from , y_to = y - y_height / zoom , y + y_height / zoom
62
+ x_width , y_height = 1.5 , 1.5 * height / width
63
+ x_from , x_to = x - x_width / zoom , x + x_width / zoom
64
+ y_from , y_to = y - y_height / zoom , y + y_height / zoom
45
65
46
66
# Here the actual algorithm starts
47
67
x = np .linspace (x_from , x_to , width ).reshape ((1 , width ))
48
68
y = np .linspace (y_from , y_to , height ).reshape ((height , 1 ))
49
- z = x + 1j * y
69
+ z = x + 1j * y
50
70
51
71
# Initialize z to all zero
52
72
c = np .full (z .shape , c )
@@ -57,16 +77,26 @@ def julia(width: int, height: int, *,
57
77
# To keep track on which points did not converge so far
58
78
m = np .full (c .shape , True , dtype = bool )
59
79
for i in range (max_iterations ):
60
- z [m ] = z [m ]** 2 + c [m ]
80
+ z [m ] = z [m ] ** 2 + c [m ]
61
81
m [np .abs (z ) > 2 ] = False
62
82
div_time [m ] = i
63
83
64
84
return div_time
65
85
86
+
66
87
Range = Tuple [float , float ]
67
88
68
- def newton (width : int , height : int , * ,
69
- p : Polynomial , a : complex , xr : Range = (- 2.5 , 1 ), yr : Range = (- 1 , 1 ), max_iterations : int = 100 ) -> (np .array , np .array ):
89
+
90
+ def newton (
91
+ width : int ,
92
+ height : int ,
93
+ * ,
94
+ p : Polynomial ,
95
+ a : complex ,
96
+ xr : Range = (- 2.5 , 1 ),
97
+ yr : Range = (- 1 , 1 ),
98
+ max_iterations : int = 100
99
+ ) -> (np .array , np .array ):
70
100
""" """
71
101
# To make navigation easier we calculate these values
72
102
x_from , x_to = xr
@@ -75,7 +105,7 @@ def newton(width: int, height: int, *,
75
105
# Here the actual algorithm starts
76
106
x = np .linspace (x_from , x_to , width ).reshape ((1 , width ))
77
107
y = np .linspace (y_from , y_to , height ).reshape ((height , 1 ))
78
- z = x + 1j * y
108
+ z = x + 1j * y
79
109
80
110
# Compute the derivative
81
111
dp = p .deriv ()
@@ -97,10 +127,12 @@ def newton(width: int, height: int, *,
97
127
r = np .full (a .shape , 0 , dtype = int )
98
128
99
129
for i in range (max_iterations ):
100
- z [m ] = z [m ] - a [m ]* p (z [m ])/ dp (z [m ])
130
+ z [m ] = z [m ] - a [m ] * p (z [m ]) / dp (z [m ])
101
131
102
132
for j , root in enumerate (roots ):
103
- converged = (np .abs (z .real - root .real ) < epsilon ) & (np .abs (z .imag - root .imag ) < epsilon )
133
+ converged = (np .abs (z .real - root .real ) < epsilon ) & (
134
+ np .abs (z .imag - root .imag ) < epsilon
135
+ )
104
136
m [converged ] = False
105
137
r [converged ] = j + 1
106
138
0 commit comments