Add a Rössler attractor

scene_diverge.py

@@ -11,6 +11,13 @@ def lorenz_system(t, state, sigma=10, rho=28, beta=8 / 3):
     dzdt = x * y - beta * z
     return [dxdt, dydt, dzdt]
 
+def roessler_system(t, state, v=3.0, a=0.2, b=0.2, c=5.7):
+    x, y, z = state
+    dxdt = v * -(y + z)
+    dydt = v * (x + a * y)
+    dzdt = v * (b + z * (x - c))
+    return [dxdt, dydt, dzdt]
+
 
 def ode_solution_points(function, state0, time, dt=0.01):
     solution = solve_ivp(
@@ -68,3 +75,53 @@ class LorenzAttractor(ThreeDScene):
             Create(curve, run_time=run_time, rate_func=linear)
             for curve in curves
         ))
+
+class RoesslerAttractor(ThreeDScene):
+    def construct(self):
+        revolutions = 1.0
+        camera_rotation_rate = 0.2
+        warm_up_time = 0.2
+        dt = 0.002
+        epsilon = 0.0025
+        states = [[1.0 + n * epsilon, 0.0, 0.0] for n in range(2)]
+        colors = color_gradient([BLUE, GREEN], len(states))
+
+        self.set_camera_orientation(phi=5 * PI / 12, theta=-6 * PI / 12)
+        self.begin_ambient_camera_rotation(rate=camera_rotation_rate)
+        axes = ThreeDAxes(
+            x_range=(-30, 30, 5),
+            y_range=(-30, 30, 5),
+            z_range=(-15, 35, 5),
+            x_length=16,
+            y_length=16,
+            z_length=12,
+        )
+        axes.center()
+        self.add(axes)
+
+        run_time = 2.0 * pi * revolutions / camera_rotation_rate
+        warm_up_len = int(warm_up_time / dt)
+
+        curves = VGroup()
+        for state, color in zip(states, colors):
+            points = ode_solution_points(
+                roessler_system, state, time=warm_up_time + run_time, dt=dt
+            )
+            curve = VMobject().set_points_as_corners(axes.c2p(points[warm_up_len:]))
+            curve.set_stroke(color)
+            curves.add(curve)
+
+        dots = Group(*(Dot3D(color=color, radius=0.1) for color in colors))
+
+        globals().update(locals())
+        def update_dots(dots):
+            for dot, curve in zip(dots, curves):
+                dot.move_to(curve.get_end())
+
+        dots.add_updater(update_dots)
+
+        self.add(dots)
+        self.play(*(
+            Create(curve, run_time=run_time, rate_func=linear)
+            for curve in curves
+        ))