我正在尝试通过使用 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 numerical 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 笔记本中使用此代码，当我执行最后一部分时：
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
我确实使用了 float32 但没有帮助。我的项目要求我对双摆进行模拟，然后使用坐标来创建概率密度，因此拥有很多点确实对模拟有很大帮助。
我想要尽可能多的点，但无法弄清楚如何做到这一点。 任何形式的帮助都会非常好:)-
64位系统也许能够在numpy.memmap中映射数组
但是如果系统是32位的话最大内存文件达到2GiB所以29.8GiB