iframe-proxy | Sunbelt Computer Software

History

220 lines (188 loc) · 6.44 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

"""Author Alexandre De Zotti

Draws Julia sets of quadratic polynomials and exponential maps.

More specifically, this iterates the function a fixed number of times

then plots whether the absolute value of the last iterate is greater than

a fixed threshold (named "escape radius"). For the exponential map this is not

really an escape radius but rather a convenient way to approximate the Julia

set with bounded orbits.

The examples presented here are:

- The Cauliflower Julia set, see e.g.

https://en.wikipedia.org/wiki/File:Julia_z2%2B0,25.png

- Other examples from https://en.wikipedia.org/wiki/Julia_set

- An exponential map Julia set, ambiantly homeomorphic to the examples in

https://www.math.univ-toulouse.fr/~cheritat/GalII/galery.html

and

https://ddd.uab.cat/pub/pubmat/02141493v43n1/02141493v43n1p27.pdf

Remark: Some overflow runtime warnings are suppressed. This is because of the

way the iteration loop is implemented, using numpy's efficient computations.

Overflows and infinites are replaced after each step by a large number.

"""

import warnings

from collections.abc import Callable

from typing import Any

import numpy

from matplotlib import pyplot

c_cauliflower = 0.25 + 0.0j

c_polynomial_1 = -0.4 + 0.6j

c_polynomial_2 = -0.1 + 0.651j

c_exponential = -2.0

nb_iterations = 56

window_size = 2.0

nb_pixels = 666

def eval_exponential(c_parameter: complex, z_values: numpy.ndarray) -> numpy.ndarray:

"""

Evaluate $e^z + c$.

>>> eval_exponential(0, 0)

1.0

>>> abs(eval_exponential(1, numpy.pi*1.j)) < 1e-15

True

>>> abs(eval_exponential(1.j, 0)-1-1.j) < 1e-15

True

"""

return numpy.exp(z_values) + c_parameter

def eval_quadratic_polynomial(

c_parameter: complex, z_values: numpy.ndarray

) -> numpy.ndarray:

"""

>>> eval_quadratic_polynomial(0, 2)

>>> eval_quadratic_polynomial(-1, 1)

>>> round(eval_quadratic_polynomial(1.j, 0).imag)

>>> round(eval_quadratic_polynomial(1.j, 0).real)

"""

return z_values * z_values + c_parameter

def prepare_grid(window_size: float, nb_pixels: int) -> numpy.ndarray:

"""

Create a grid of complex values of size nb_pixels*nb_pixels with real and

imaginary parts ranging from -window_size to window_size (inclusive).

Returns a numpy array.

>>> prepare_grid(1,3)

array([[-1.-1.j, -1.+0.j, -1.+1.j],

[ 0.-1.j, 0.+0.j, 0.+1.j],

[ 1.-1.j, 1.+0.j, 1.+1.j]])

"""

x = numpy.linspace(-window_size, window_size, nb_pixels)

x = x.reshape((nb_pixels, 1))

y = numpy.linspace(-window_size, window_size, nb_pixels)

y = y.reshape((1, nb_pixels))

return x + 1.0j * y

def iterate_function(

eval_function: Callable[[Any, numpy.ndarray], numpy.ndarray],

function_params: Any,

nb_iterations: int,

z_0: numpy.ndarray,

infinity: float | None = None,

) -> numpy.ndarray:

"""

Iterate the function "eval_function" exactly nb_iterations times.

The first argument of the function is a parameter which is contained in

function_params. The variable z_0 is an array that contains the initial

values to iterate from.

This function returns the final iterates.

>>> iterate_function(eval_quadratic_polynomial, 0, 3, numpy.array([0,1,2])).shape

(3,)

>>> numpy.round(iterate_function(eval_quadratic_polynomial,

... 0,

... 3,

... numpy.array([0,1,2]))[0])

>>> numpy.round(iterate_function(eval_quadratic_polynomial,

... 0,

... 3,

... numpy.array([0,1,2]))[1])

(1+0j)

>>> numpy.round(iterate_function(eval_quadratic_polynomial,

... 0,

... 3,

... numpy.array([0,1,2]))[2])

(256+0j)

"""

z_n = z_0.astype("complex64")

for _ in range(nb_iterations):

z_n = eval_function(function_params, z_n)

if infinity is not None:

numpy.nan_to_num(z_n, copy=False, nan=infinity)

z_n[abs(z_n) == numpy.inf] = infinity

return z_n

def show_results(

function_label: str,

function_params: Any,

escape_radius: float,

z_final: numpy.ndarray,

) -> None:

"""

Plots of whether the absolute value of z_final is greater than

the value of escape_radius. Adds the function_label and function_params to

the title.

>>> show_results('80', 0, 1, numpy.array([[0,1,.5],[.4,2,1.1],[.2,1,1.3]]))

"""

abs_z_final = (abs(z_final)).transpose()

abs_z_final[:, :] = abs_z_final[::-1, :]

pyplot.matshow(abs_z_final < escape_radius)

pyplot.title(f"Julia set of ${function_label}$, $c={function_params}$")

pyplot.show()

def ignore_overflow_warnings() -> None:

"""

Ignore some overflow and invalid value warnings.

>>> ignore_overflow_warnings()

"""

warnings.filterwarnings(

"ignore", category=RuntimeWarning, message="overflow encountered in multiply"

)

warnings.filterwarnings(

"ignore",

category=RuntimeWarning,

message="invalid value encountered in multiply",

)

warnings.filterwarnings(

"ignore", category=RuntimeWarning, message="overflow encountered in absolute"

)

warnings.filterwarnings(

"ignore", category=RuntimeWarning, message="overflow encountered in exp"

)

if __name__ == "__main__":

z_0 = prepare_grid(window_size, nb_pixels)

ignore_overflow_warnings() # See file header for explanations

nb_iterations = 24

escape_radius = 2 * abs(c_cauliflower) + 1

z_final = iterate_function(

eval_quadratic_polynomial,

c_cauliflower,

nb_iterations,

z_0,

infinity=1.1 * escape_radius,

)

show_results("z^2+c", c_cauliflower, escape_radius, z_final)

nb_iterations = 64

escape_radius = 2 * abs(c_polynomial_1) + 1

z_final = iterate_function(

eval_quadratic_polynomial,

c_polynomial_1,

nb_iterations,

z_0,

infinity=1.1 * escape_radius,

)

show_results("z^2+c", c_polynomial_1, escape_radius, z_final)

nb_iterations = 161

escape_radius = 2 * abs(c_polynomial_2) + 1

z_final = iterate_function(

eval_quadratic_polynomial,

c_polynomial_2,

nb_iterations,

z_0,

infinity=1.1 * escape_radius,

)

show_results("z^2+c", c_polynomial_2, escape_radius, z_final)

nb_iterations = 12

escape_radius = 10000.0

z_final = iterate_function(

eval_exponential,

c_exponential,

nb_iterations,

z_0 + 2,

infinity=1.0e10,

)

show_results("e^z+c", c_exponential, escape_radius, z_final)

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

julia_sets.py

Latest commit

History

julia_sets.py

File metadata and controls