Python / Scientific and math computing
Exploring Arithmetic of Dynamical Systems with Python
Contents
Introduction to Dynamical Systems
Dynamical systems are mathematical models that describe how a system changes over time. In the context of arithmetic, we focus on discrete dynamical systems where the state evolves in discrete time steps.
def iterate(f, x0, n):
"""Iterate function f n times starting from x0."""
x = x0
for _ in range(n):
x = f(x)
return x
# Example: f(x) = x^2 - 1
f = lambda x: x**2 - 1
x0 = 0.5
print(iterate(f, x0, 10))
Fixed Points
Fixed points are values that remain unchanged when the function is applied. They are crucial in understanding the long-term behavior of dynamical systems.
import numpy as np
def find_fixed_points(f, x_range):
x = np.linspace(x_range[0], x_range[1], 1000)
y = f(x)
fixed_points = x[np.abs(x - y) < 1e-6]
return fixed_points
# Example: f(x) = x^2 - 1
f = lambda x: x**2 - 1
fixed_points = find_fixed_points(f, (-2, 2))
print("Fixed points:", fixed_points)
Orbits
An orbit is the sequence of points generated by repeatedly applying the function to an initial value.
def plot_orbit(f, x0, n):
import matplotlib.pyplot as plt
orbit = [x0]
for _ in range(n):
x0 = f(x0)
orbit.append(x0)
plt.plot(range(n+1), orbit, 'b.-')
plt.xlabel('Iteration')
plt.ylabel('Value')
plt.title('Orbit of x0')
plt.show()
# Example: f(x) = 3.2 * x * (1 - x)
f = lambda x: 3.2 * x * (1 - x)
plot_orbit(f, 0.4, 50)
Periodic Points
Periodic points are values that return to themselves after a finite number of iterations.
def find_periodic_points(f, x_range, period):
x = np.linspace(x_range[0], x_range[1], 1000)
y = x
for _ in range(period):
y = f(y)
periodic_points = x[np.abs(x - y) < 1e-6]
return periodic_points
# Example: f(x) = 4x(1-x)
f = lambda x: 4 * x * (1 - x)
periodic_points = find_periodic_points(f, (0, 1), 2)
print("Period-2 points:", periodic_points)
Bifurcation Diagrams
Bifurcation diagrams show how the long-term behavior of a system changes as a parameter varies.
def plot_bifurcation(f, r_range, x0, n, n_discard):
import matplotlib.pyplot as plt
r_values = np.linspace(r_range[0], r_range[1], 1000)
x_values = []
for r in r_values:
x = x0
for _ in range(n_discard):
x = f(x, r)
for _ in range(n):
x = f(x, r)
x_values.append((r, x))
r_plot, x_plot = zip(*x_values)
plt.plot(r_plot, x_plot, ',k', alpha=0.1)
plt.xlabel('r')
plt.ylabel('x')
plt.title('Bifurcation Diagram')
plt.show()
# Example: Logistic Map
f = lambda x, r: r * x * (1 - x)
plot_bifurcation(f, (2.5, 4), 0.5, 100, 100)
Lyapunov Exponents
Lyapunov exponents measure the rate of separation of infinitesimally close trajectories, indicating chaos.
def lyapunov_exponent(f, x0, n, epsilon=1e-6):
x = x0
lyap = 0
for _ in range(n):
x_perturbed = x + epsilon
lyap += np.log(abs((f(x_perturbed) - f(x)) / epsilon))
x = f(x)
return lyap / n
# Example: Logistic Map
f = lambda x: 3.9 * x * (1 - x)
lyap = lyapunov_exponent(f, 0.5, 1000)
print(f"Lyapunov exponent: {lyap:.6f}")
Fractals and Self-Similarity
Some dynamical systems exhibit fractal structures and self-similarity.
def mandelbrot(h, w, max_iter):
y, x = np.ogrid[-1.4:1.4:h*1j, -2:0.8:w*1j]
c = x + y*1j
z = c
divtime = max_iter + np.zeros(z.shape, dtype=int)
for i in range(max_iter):
z = z**2 + c
diverge = z*np.conj(z) > 2**2
div_now = diverge & (divtime == max_iter)
divtime[div_now] = i
z[diverge] = 2
return divtime
import matplotlib.pyplot as plt
plt.imshow(mandelbrot(500, 750, 50), cmap='hot', extent=[-2, 0.8, -1.4, 1.4])
plt.title('Mandelbrot Set')
plt.show()
Symbolic Dynamics
Symbolic dynamics represents the behavior of a dynamical system using sequences of symbols.
def itinerary(x0, n, partition_point):
x = x0
itinerary = []
for _ in range(n):
if x < partition_point:
itinerary.append('0')
else:
itinerary.append('1')
x = 4 * x * (1 - x) # Logistic map with r = 4
return ''.join(itinerary)
# Example
x0 = 0.7
partition_point = 0.5
print(itinerary(x0, 20, partition_point))
Recurrence Plots
Recurrence plots visualize the recurrence of states in phase space.
def recurrence_plot(x, epsilon):
dist = np.abs(x[:, None] - x)
return (dist < epsilon).astype(int)
# Generate time series
t = np.linspace(0, 20, 1000)
x = np.sin(t) + 0.1 * np.random.randn(1000)
rp = recurrence_plot(x, 0.2)
plt.imshow(rp, cmap='binary')
plt.title('Recurrence Plot')
plt.show()
Poincaré Maps
Poincaré maps reduce the dimensionality of a continuous dynamical system by examining its intersection with a lower-dimensional subspace.
from scipy.integrate import odeint
def rossler_system(state, t, a, b, c):
x, y, z = state
return [-y - z, x + a*y, b + z*(x - c)]
def poincare_map(trajectory, plane='y'):
x, y, z = trajectory.T
if plane == 'y':
return x[np.where(np.diff(np.sign(y)) > 0)[0]], z[np.where(np.diff(np.sign(y)) > 0)[0]]
a, b, c = 0.2, 0.2, 5.7
t = np.linspace(0, 500, 50000)
trajectory = odeint(rossler_system, [1, 1, 1], t, args=(a, b, c))
x_poincare, z_poincare = poincare_map(trajectory)
plt.plot(x_poincare, z_poincare, ',k')
plt.title('Poincaré Map of Rössler System')
plt.xlabel('x')
plt.ylabel('z')
plt.show()
Basin of Attraction
The basin of attraction is the set of initial conditions that evolve towards a particular attractor.
def basin_of_attraction(f, x_range, y_range, n_points, n_iter):
x = np.linspace(x_range[0], x_range[1], n_points)
y = np.linspace(y_range[0], y_range[1], n_points)
X, Y = np.meshgrid(x, y)
Z = X + 1j*Y
for _ in range(n_iter):
Z = f(Z)
return abs(Z) < 2
# Example: z^2 + c, where c = -0.4 + 0.6j
c = -0.4 + 0.6j
f = lambda z: z**2 + c
basin = basin_of_attraction(f, (-2, 2), (-2, 2), 500, 50)
plt.imshow(basin, extent=[-2, 2, -2, 2], cmap='binary')
plt.title('Basin of Attraction')
plt.xlabel('Re(z)')
plt.ylabel('Im(z)')
plt.show()
Stability Analysis
Stability analysis determines how small perturbations affect the long-term behavior of a system.
def stability_analysis(f, fixed_point, epsilon=1e-6):
derivative = (f(fixed_point + epsilon) - f(fixed_point)) / epsilon
if abs(derivative) < 1:
return "Stable"
elif abs(derivative) > 1:
return "Unstable"
else:
return "Neutral"
# Example: f(x) = x^2 - 1
f = lambda x: x**2 - 1
fixed_points = [1.618033988749895, -0.618033988749895] # Golden ratio and its negative
for fp in fixed_points:
stability = stability_analysis(f, fp)
print(f"Fixed point {fp:.6f} is {stability}")
Chaos and Sensitivity to Initial Conditions
Chaotic systems exhibit sensitive dependence on initial conditions, where small changes lead to drastically different outcomes.
def plot_sensitivity(f, x1, x2, n):
trajectory1 = [x1]
trajectory2 = [x2]
for _ in range(n):
x1 = f(x1)
x2 = f(x2)
trajectory1.append(x1)
trajectory2.append(x2)
plt.plot(range(n+1), trajectory1, 'b-', label='x0 = %.6f' % trajectory1[0])
plt.plot(range(n+1), trajectory2, 'r-', label='x0 = %.6f' % trajectory2[0])
plt.xlabel('Iteration')
plt.ylabel('Value')
plt.title('Sensitivity to Initial Conditions')
plt.legend()
plt.show()
# Example: Logistic Map in chaotic regime
f = lambda x: 3.9 * x * (1 - x)
plot_sensitivity(f, 0.5, 0.500001, 50)
Additional Resources
For further exploration of dynamical systems and their arithmetic properties, consider the following resources:
- "Chaos: An Introduction to Dynamical Systems" by Kathleen T. Alligood, Tim D. Sauer, and James A. Yorke ArXiv: https://arxiv.org/abs/nlin/0001016
- "Introduction to the Modern Theory of Dynamical Systems" by Anatole Katok and Boris Hasselblatt ArXiv: https://arxiv.org/abs/math/9903105
- "Ergodic Theory and Dynamical Systems" journal https://www.cambridge.org/core/journals/ergodic-theory-and-dynamical-systems
These resources provide in-depth discussions on various aspects of dynamical systems, including their arithmetic properties and applications in different fields of mathematics and science.