我试图通过使用numpy在python中求解必要的微分方程来模拟一个双摆系统。
下面是我的代码：

import numpy as np
import sympy as smp
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import PillowWriter
from sympy.physics.mechanics import dynamicsymbols, init_vprinting
import pandas as pd
import matplotlib.pyplot as plt
t, g = smp.symbols('t g')
m1, m2 = smp.symbols('m1 m2')
L1, L2 = smp.symbols('L1, L2')
the1, the2 = smp.symbols(r'\theta_1, \theta_2', cls=smp.Function)
the1 = the1(t)
the2 = the2(t)
the1_d = the1.diff(t)
the2_d = the2.diff(t)
the1_dd = the1_d.diff(t)
the2_dd = the2_d.diff(t)
x1 = L1*smp.sin(the1)
y1 = -L1*smp.cos(the1)
x2 = L1*smp.sin(the1)+L2*smp.sin(the2)
y2 = -L1*smp.cos(the1)-L2*smp.cos(the2)
# Kinetic
T1 = 1/2 * m1 * (smp.diff(x1, t)**2 + smp.diff(y1, t)**2)
T2 = 1/2 * m2 * (smp.diff(x2, t)**2 + smp.diff(y2, t)**2)
T = T1+T2
# Potential
V1 = m1*g*y1
V2 = m2*g*y2
V = V1 + V2
# Lagrangian
L = T-V
#using the  the euler-langragian
LE1 = smp.diff(L, the1) - smp.diff(smp.diff(L, the1_d), t).simplify()
LE2 = smp.diff(L, the2) - smp.diff(smp.diff(L, the2_d), t).simplify()

#now since we have two linear equations with two unknowns that is the second order s´dervative we can slve the equations 
sols = smp.solve([LE1, LE2], (the1_dd, the2_dd),
                simplify=False, rational=False)

#we have a symbolic expression above which can be converted into a numberical function provided that the values are given and the differential equations are defined
dz1dt_f = smp.lambdify((t,g,m1,m2,L1,L2,the1,the2,the1_d,the2_d), sols[the1_dd])
dz2dt_f = smp.lambdify((t,g,m1,m2,L1,L2,the1,the2,the1_d,the2_d), sols[the2_dd])
dthe1dt_f = smp.lambdify(the1_d, the1_d)
dthe2dt_f = smp.lambdify(the2_d, the2_d)

def dSdt(S, t, g, m1, m2, L1, L2):
    the1, z1, the2, z2 = S
    return [
        dthe1dt_f(z1),
        dz1dt_f(t, g, m1, m2, L1, L2, the1, the2, z1, z2),
        dthe2dt_f(z2),
        dz2dt_f(t, g, m1, m2, L1, L2, the1, the2, z1, z2),
    ]

t = np.linspace(0, 40000000,1000000001,dtype="float32" )
g = 9.81
m1 = 50
m2=25
L1 = 1
L2=1
ans = odeint(dSdt, y0=[1.57, 0, 1.57, 0], t=t, args=(g,m1,m2,L1,L2))

现在我在jupyter notebook中使用这段代码，当我执行最后一部分时：

t = np.linspace(0, 40000000,1000000001,dtype="float32" )
g = 9.81
m1 = 50
m2=25
L1 = 1
L2=1
ans = odeint(dSdt, y0=[1.57, 0, 1.57, 0], t=t, args=(g,m1,m2,L1,L2))

我得到这个错误：

MemoryError                               Traceback (most recent call last)
Cell In [19], line 7
      5 L1 = 1
      6 L2=1
----> 7 ans = odeint(dSdt, y0=[1.57, 0, 1.57, 0], t=t, args=(g,m1,m2,L1,L2))

File d:\python\lib\site-packages\scipy\integrate\_odepack_py.py:241, in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    239 t = copy(t)
    240 y0 = copy(y0)
--> 241 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
    242                          full_output, rtol, atol, tcrit, h0, hmax, hmin,
    243                          ixpr, mxstep, mxhnil, mxordn, mxords,
    244                          int(bool(tfirst)))
    245 if output[-1] < 0:
    246     warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

MemoryError: Unable to allocate 29.8 GiB for an array with shape (1000000001, 4) and data type float64

我确实使用了Float 32，但它没有帮助。我的项目要求我做双摆的模拟，然后使用坐标创建概率密度，所以有很多点确实对模拟有很大帮助。
我想有尽可能多的点，但不能想出如何做到这一点。任何形式的帮助都会非常好：）-