Анимация двойного маятника (питон)

Question

Я работаю над приближением движения системы с двойным маятником, используя несколько решателей ODE, и мне очень хотелось бы иметь возможность показывать 4 разные анимации одновременно, используя вспомогательные сюжеты. Я начинающий программист и не могу понять, как использовать вспомогательные сюжеты для одновременного отображения анимации.

Обратите внимание, что мои решатели ODE находятся в другом файле py. Спасибо!

from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from solvers import euler, trapezoid

# Constant Values
g = 9.8   # acceleration due to gravity, in m/s^2
L1 = 1.0  # length of pendulum 1 in m
L2 = 1.0  # length of pendulum 2 in m
m1 = 1.0  # mass of pendulum 1 in kg
m2 = 1.0  # mass of pendulum 2 in kg

# Initial Conditions
  # th1 and th2 are the angles (degrees)
  # w1 and w2 are the angular velocities (degrees per second)
th1 = 0
w1 = 0.5
th2 = 0
w2 = 0.5

# Convert IC's into radians and radians per second
cond = np.radians([th1, w1, th2, w2]) # the current state

# Define the f(theta,omega) functions for each derivative
def f(cond):
    y=np.zeros([len(cond)])

    # Store values that repeat often
    COS = cos(cond[0]-cond[2])
    SIN = sin(cond[0]-cond[2])

    # d(theta_1)/dt = w1
    y[0] = cond[1]

    # d(w_1)/dt
    y[1] = (m2 * g * sin(cond[2]) * COS
            - m2 * SIN * (L1 * cond[1] ** 2 * COS + L2 * cond[3] ** 2)
            - (m1 + m2) * g * sin(cond[0])) \
           / (L1 * (m1 + m2 * SIN ** 2))

    # d(theta_2)/dt = w2
    y[2] = cond[3]

    # d(w_2)/dt
    y[3] = ((m1 + m2)*(L1 * cond[1]**2 * SIN - g * sin(cond[2])
            + g*sin(cond[0]) * COS) + m2 * L2 * cond[3]**2 * SIN * COS) \
            / (L2*(m1 + m2 * SIN**2))

    return y #return the f(theta,omega)


# Convert the solution into cartesian coordinates
def cartesian(sol):
    xy = np.zeros((len(sol),4))
    xy[:,0] = L1 * sin(sol[:, 0])     # x1
    xy[:,1] = -L1 * cos(sol[:, 0])    # y1

    xy[:,2] = L2 * sin(sol[:, 2]) + xy[:,0]
    xy[:,3] = -L2 * cos(sol[:, 2]) + xy[:,1]

    return xy

# Animation
def animate(y,h, name):
    x1 = y[:,0]
    y1 = y[:,1]
    x2 = y[:,2]
    y2 = y[:,3]


    fig = plt.figure()
    ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
    ax.grid()

    line, = ax.plot([], [], 'o-', lw=2)
    time_template = 'time = %.1fs'
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

    def init():
        line.set_data([], [])
        time_text.set_text('')
        return line, time_text

    def animate(i):
        thisx = [0, x1[i], x2[i]]
        thisy = [0, y1[i], y2[i]]

        line.set_data(thisx, thisy)
        time_text.set_text(time_template % (i * h))
        return line, time_text

    ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y)),
                              interval=30, blit=True, init_func=init)

    plt.title(name)

    return plt.show()


# Create a time array (20 seconds, n = 400 => h = .05)
tmax = 20
n = 400
t = np.linspace(0, tmax, n)
h = (t[len(t)-1]-t[0])/len(t)


# Store the Approximations
eul = cartesian(euler(f, cond, t))
trap = cartesian(trapezoid(f,cond, t))


# Show Animation
animate(eul,h, 'Euler\'s Method')
animate(trap, h, 'Trapezoid Method')