アニメーションの実装方法について教えてください。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

N = 2
K = 5
w1 = np.random.standard_cauchy()
w2 = np.random.standard_cauchy()
T = 50
n = 500
h = T / n
ts = np.arange(0, T, h)

f1 = lambda theta1,theta2,t : w1-K*np.sin(theta1-theta2)/N
f2 = lambda theta1,theta2,t : w2-K*np.sin(theta2-theta1)/N

theta1 = np.pi / 2
theta2 = 3 * np.pi / 2

fig, ax = plt.subplots()
ax.set_aspect("equal")

ims = []
for t in ts:
	# 描画
	x = np.cos([theta1, theta2])
	y = np.sin([theta1, theta2])

	im = ax.scatter(x, y, facecolor="None", edgecolors="r")
	ims.append([im]) # 注意：リストで渡す必要があります。

	# 更新
	k1 = h * f1(theta1, theta2, t)
	l1 = h * f2(theta1, theta2, t)
	k2 = h * f1(theta1 + k1 / 2, theta2 + l1 / 2, t + h / 2)
	l2 = h * f2(theta1 + k1 / 2, theta2 + l1 / 2, t + h / 2)
	k3 = h * f1(theta1 + k2 / 2, theta2 + l2 / 2, t + h / 2)
	l3 = h * f2(theta1 + k2 / 2, theta2 + l2 / 2, t + h / 2)
	k4 = h * f1(theta1 + k3, theta2 + l3, t + h)
	l4 = h * f2(theta1 + k3, theta2 + l3, t + h)

	theta1 += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4 )
	theta2 += 1/6 * (l1 + 2 * l2 + 2 * l3 + l4 )

anim = animation.ArtistAnimation(fig, ims, interval=100)
plt.show()

import math
import random
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

# コーシー分布
def rand_cauchy(mu, gamma) :
	arg = math.pi * random.uniform(-0.5, 0.5)
	return mu + gamma * math.tan(arg)

# 微分方程式の右辺
def func(t, theta, **params):
	K = params.get("K", 0.0)
	omega = params.get("omega", np.zeros_like(theta))
	Ti, Tj = np.meshgrid(theta, theta)
	return omega + K * np.average(np.sin(Tj - Ti), axis=0)

# 4次ルンゲクッタ
def rk4(dt, t, Xs, func, **params):
	k1 = dt * func(t, Xs, **params)
	k2 = dt * func(t + dt / 2, Xs + k1 / 2, **params)
	k3 = dt * func(t + dt / 2, Xs + k2 / 2, **params)
	k4 = dt * func(t + dt, Xs + k3, **params)
	# 更新
	new_t = t + dt
	newXs = Xs + (k1 + 2 * k2 + 2 * k3 + k4) / 6
	return new_t, newXs

# パラメータ
N = 100
K = 10
T = 5
dt = 0.01
ts = np.arange(0, T + dt, dt)

# 初期値
# omega = np.array([np.random.standard_cauchy() for _ in range(N)])
omega = np.array([rand_cauchy(4, 1) for _ in range(N)])
theta = np.array([np.random.uniform(0.0, 2 * np.pi) for _ in range(N)])

# 計算の実行
fig, ax = plt.subplots()
ax.set_aspect("equal")

ims = []
for t in ts:
	# 描画
	x = np.cos(theta)
	y = np.sin(theta)

	im = ax.scatter(x, y, facecolor="None", edgecolors="r")
	ims.append([im]) # 注意：リストで渡す必要があります。

	# 更新
	_, theta = rk4(dt, t, theta, func, K=K, omega=omega)

ani = animation.ArtistAnimation(fig, ims, interval=100)
# ani.save("./animation.gif")
plt.show()