In [2]: 1 import time
2 import numpy as np
3 import sympy as sy
4 import [Link] as plt
5 from [Link] import interp1d
6 import [Link] as cm
7 import [Link] as mcolors
8 from collections import defaultdict
9 import warnings
10 [Link]('ignore')
11 [Link]['[Link]'] = True
12 [Link]['[Link]'] = 20
13 [Link]['[Link]'] = 30
14 from [Link] import odeint
𝑥¨ + 𝛿𝑥˙ − 𝑥 + 𝑥3 = 𝐹 cos 𝜔𝑡
In [3]: 1 def double_well_oscillator( y, t, delta, omega, F):
2 x, v = y
3 dxdt = v
4 dvdt = - delta * v + x - x**3 + F * [Link](omega * t)
5 return [dxdt, dvdt]
(Basins for the unforced oscillator)
Sketch the basins for the weakly damped double-well oscillator (12.5.1) in the unforced case when 𝐹 = 0.
How does their shape depend on the size of the damping? What happens to the basins as the damping tends to zero?
What implications does this have for the predictability of the unforced system?
In [10]: 1 F = 0
2 omega = 1
3 T = 2 * [Link] / omega
4 t = [Link](0, 73*T, 1001) # Time vector
5 delta_values = [0.25, 0.125, 0.0625] # Damping values to explore
6
7 x_vals = [Link](-2.5, 2.5, 900) # Initial positions
8 v_vals = [Link](-2.5, 2.5, 900) # Initial velocities
9 X, V = [Link]( x_vals, v_vals)
10 initial_conditions = np.column_stack([[Link](), [Link]()])
11 start_time = [Link]()
12 fig, axs = [Link](1, len(delta_values), figsize=(15, 5), sharex=True, sharey=True)
13 for idx, delta in enumerate(delta_values):
14 basins = np.zeros_like(X)
15 for i, y0 in enumerate(initial_conditions):
16 sol = odeint(double_well_oscillator, y0, t, args=(delta, omega, F))
17 x_final = sol[-1, 0] # Final x position after integration
18 if x_final < 0:
19 [Link]()[i] = -1 # Left well
20 else:
21 [Link]()[i] = 1 # Right well
22
23 ax = axs[idx]
24 [Link](X, V, basins, level=[-2.5,0,2.5], cmap='jet')
25 ax.set_title(f'Damping: $\delta = {delta}$')
26 ax.set_xlabel('x (Position)')
27 ax.set_ylabel('v (Velocity)')
28 [Link](r'Basins of Attraction for the Double-Well Oscillator , Parameter : $F=0,\omega=1$')
29 plt.tight_layout()
30 [Link]()
31
32 end_time = [Link]()
33
34 elapsed_time = end_time - start_time # Calculate the total elapsed time
35
36 # Convert the elapsed time into hours, minutes, and seconds
37 hours = int(elapsed_time // 3600)
38 minutes = int((elapsed_time % 3600) // 60)
39 seconds = elapsed_time % 60
40 print(f"Total execution time: {hours} hours, {minutes} minutes, and {seconds:.4f} seconds")
Total execution time: 12 hours, 42 minutes, and 26.4108 seconds
The basins become thinner as the damping decreases, and they appear to wrap around each other more and more as 𝛿 → 0 (tends to zero).
This will make for sensitive dependence on initial conditions, which means that accurate prediction of which basin an initial condition will settle into is virtually impossible unless the
initial condition is practically settled into one of the basins already.
(Basins for the forced oscillator)
Sketch the basins for the weakly damped double-well oscillator (12.5.1) in the forced case when 𝐹 = 0.25, 𝜔 = 1.
In [5]: 1 F = 0.25
2 delta = 0.25
3 omega = 1
4 T = 2 * [Link] / omega
5 t = [Link](0, 73*T, 1000) # Time vector
6 x_vals = [Link](-2.5, 2.5, 900) # Initial positions
7 v_vals = [Link](-2.5, 2.5, 900) # Initial velocities
8 X, V = [Link]( x_vals, v_vals)
9 initial_conditions = np.column_stack([[Link](), [Link]()])
10
11 fig, ax = [Link]( figsize=(30, 27), sharex=True, sharey=True)
12 start_time = [Link]()
13 basins = np.zeros_like(X)
14
15 for i, y0 in enumerate(initial_conditions):
16
17 sol = odeint(double_well_oscillator, y0, t, args=(delta, omega, F))
18 x_final = sol[-1, 0] # Final x position after integration
19 if x_final < 0:
20 [Link]()[i] = -1 # Left well
21 else:
22 [Link]()[i] = 1 # Right well
23
24 [Link](X, V, basins, levels=[-2.5, 0, 2.5], cmap='jet')
25 ax.set_title(f'Damping: $\delta = {delta}$')
26 ax.set_xlabel('x (Position)')
27 ax.set_ylabel('v (Velocity)')
28 [Link]('Basins of Attraction for the Double-Well Oscillator')
29 plt.tight_layout()
30 [Link]()
31
32 end_time = [Link]()
33
34 elapsed_time = end_time - start_time # Calculate the total elapsed time
35
36 # Convert the elapsed time into hours, minutes, and seconds
37 hours = int(elapsed_time // 3600)
38 minutes = int((elapsed_time % 3600) // 60)
39 seconds = elapsed_time % 60
40 print(f"Total execution time: {hours} hours, {minutes} minutes, and {seconds:.4f} seconds")
Total execution time: 10 hours, 15 minutes, and 44.0163 seconds
(Coexisting chaos and limit cycle)
Consider the double-well oscillator (12.5.1) with parameters 𝛿 = 0.15, 𝐹 = 0.3, and 𝜔 = 1.
Show numerically that the system has at least two coexisting attractors : a large limit cycle and a smaller strange attractor. Plot both in a Poincaré section.
In [5]: 1 delta = 0.15
2 F = 0.3
3 omega = 1.0
4 T = 2 * [Link] / omega # Period of the driving force
5 t_span = [Link](0, 5000 * T, 500000) # Integrate over 500 periods
6 poincare_times = [Link](0, 5000 * T, T) # Poincaré section times
7
8 initial_conditions = [ [0.5, 0], # Initial condition 1 (expecting large limit cycle)
9 [1.5, 0], ] # Initial condition 2 (expecting strange attractor)
10
11 solutions = [odeint(double_well_oscillator, y0, t_span, args=(delta, omega, F)) for y0 in initial_conditions]
12
13 fig, ax = [Link](figsize=(8, 6))
14 colors = ['blue', 'red']
15 labels = ['Large Limit Cycle', 'Strange Attractor']
16
17 for i, solution in enumerate(solutions):
18 x, v = [Link](poincare_times, t_span, solution[:, 0]), [Link](poincare_times, t_span, solution[:, 1])
19 [Link](x, v, s=0.5, color=colors[i], label=labels[i])
20 ax.set_title('Poincaré Section: Coexisting Attractors in Double-Well Oscillator')
21 ax.set_xlabel('x (Position)')
22 ax.set_ylabel('v (Velocity)')
23 [Link](fontsize=10)
24 [Link](True)
25 [Link]()
In [ ]: 1
