N重振り子　先端の角度のみ変えたい

# coding: utf-8
import math
import numpy as np

def create_pendulums(settings):
    pendulums = []
    for setting in settings:
        p = NthPendulum(setting["name"], setting["max_len"])
        p.set_param(setting["mass"], setting["length"])
        p.set_init_val(setting["t0"], setting["theta0"], setting["theta_dot0"])
        pendulums.append(p)
    return pendulums

class NthPendulum:
    def __init__(self, name, max_len):
        self.name = name
        self.orbit = {"max_len": max_len, "x": [], "y": []}

    def set_param(self, m, l):
        self.g = 9.8
        self.m = m
        self.l = l
        self.n = m.shape[0]

    def get_orbit(self, i):
        x, y = [], []
        for xn, yn in zip(self.orbit["x"], self.orbit["y"]):
            x.append(xn[i])
            y.append(yn[i])
        return x, y


    def set_init_val(self, t0, theta0, theta_dot0):
        self.t = t0
        self.theta = theta0
        self.theta_dot = theta_dot0

    def get_position(self):
        x =   self.l * np.sin(self.theta)
        y = - self.l * np.cos(self.theta)

        for i in range(self.n):
            if i > 0:
                x[i] += x[i-1]
                y[i] += y[i-1]

        self._save_position(x, y)

        return x, y

    def _save_position(self, x, y):
        self.orbit["x"].append(x)
        self.orbit["y"].append(y)

        if len(self.orbit["x"]) > self.orbit["max_len"]:
            self.orbit["x"].pop(0)
            self.orbit["y"].pop(0)

    def update_position(self, dt):
        new_vals = rk4(dt, self.t, self.theta, self.theta_dot, self.func)
        self.t, self.theta, self.theta_dot = new_vals

    def func(self, t, theta, theta_dot):
        I = np.ones(self.n, float)
        Z = np.zeros((self.n, self.n), float)
        A = np.matrix(Z.copy())
        B = np.matrix(Z.copy())
        E = np.matrix(I).T

        for i in range(self.n):
            for j in range(self.n):
                for k in range(max(i, j), self.n):
                    A[i, j] += self.m[k]
                    B[i, j] += self.m[k]
                if i == j:
                    A[i, j] *= self.l[j]
                    B[i, j] *= self.g * math.sin(theta[i])
                else:
                    A[i, j] *= self.l[j] * math.cos(theta[i] - theta[j])
                    B[i, j] *= self.l[j] * theta_dot[j]**2 * math.sin(theta[i] - theta[j])

        f = - np.linalg.inv(A) * B * E

        return np.array(f).flatten()

def rk4(dt, t0, x0, v0, func):
    '''
    dx/dt = v
    dv/dt = func(t, x, v)
    '''

    k1_x = v0 * dt
    k1_v = func(t0, x0, v0) * dt

    k2_x = (v0 + k1_v / 2) * dt
    k2_v = func(t0 + dt / 2, x0 + k1_x / 2, v0 + k1_v / 2) * dt

    k3_x = (v0 + k2_v / 2) * dt
    k3_v = func(t0 + dt / 2, x0 + k2_x / 2, v0 + k2_v / 2) * dt

    k4_x = (v0 + k3_v) * dt
    k4_v = func(t0 + dt, x0 + k3_x, v0 + k3_v) * dt

    t = t0 + dt
    x = x0 + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6
    v = v0 + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6

    return t, x, v

# coding: utf-8
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# <計算の設定>
t0 =  0.0 # 開始時間
te = 100.0 # 終了時間
dt = 0.05 # 時間刻み

# <アニメーションの設定>
intrvl = 40             # アニメーションの間隔
cnt = int((te - t0) / dt) # アニメーションの画像の枚数
orbit= True         # 軌道の表示の有無
max_len = 20              # 振り子の軌道のリストの長さ

# <振り子のパラメータと初期値> ここは、ファイルを読み込む形式にしてもよい

pnum = 2     # 振り子の数
mass = 1.0       # 振り子の質量
length = 1.0     # 振り子の長さ
theta_dot0 = 0.0 # 振り子の初期角速度

vnum = 2       # 初期角度のバリエーション数
a0 = np.radians(120) # 初期角度のバリエーションの範囲下限
ae = np.radians(125)# 初期角度のバリエーションの範囲上限
theta0_list=[a0,ae]

settings = [
    {
        "name":  f"$\\theta_0 = {math.degrees(theta0_list[i]):.2f}$",
        "t0": t0,
        "max_len": max_len,
        "mass": np.array([mass for j in range(pnum)]),
        "length": np.array([length for j in range(pnum)]),
        "theta0": np.array([theta0_list[i] for j in range(pnum)]),
        "theta_dot0": np.array([theta_dot0 for j in range(pnum)])
    } for i in range(len(theta0_list))
]

pendulums = create_pendulums(settings)

# <--- MAIN --->

def gen():
    # 順番に値を更新していく
    while True:
        for p in pendulums: p.update_position(dt)
        print(f"t = {p.t:.3f} sec ({p.t/te*100:.1f} %)\r", end="")
        yield pendulums


def update(pendulums):

    # 図をクリア
    ax.cla()

    t, R = 0, 0
    origin = np.array([0])

    for p in pendulums:
        t, R = p.t, max(R, 1.2 * p.l.sum())

        x, y = p.get_position()
        ox = np.concatenate([origin, x])
        oy = np.concatenate([origin, y])

        lines, = ax.plot(ox, oy)
        ax.scatter(x, y, color=lines.get_color(), label=p.name)
       # if orbit:
            # すべての振り子の軌道を表示
            #for i in range(p.n):
                #x_orbit, y_orbit = p.get_orbit(i)
                #ax.plot(x_orbit, y_orbit, color=lines.get_color(), alpha=0.5)

            # 最後の振り子の軌道のみ表示
            # x_orbit, y_orbit = p.get_orbit(p.n-1)
            # ax.plot(x_orbit, y_orbit, color=lines.get_color(), alpha=0.5)
    ax.text(
        0.05, 0.95, f"$t = {t:.3f}$ sec",
        verticalalignment='top', transform=ax.transAxes
    )

    ax.set_aspect("equal")
    ax.set_xlim(-R, R)
    ax.set_ylim(-R, R)

fig, ax = plt.subplots()
ani = FuncAnimation(fig,update,gen,interval=intrvl, save_count=cnt-1)
#ani.save("N_pendulum 3.gif")