旋转圆盘上悬挂摆的模拟
有没有人能让这段代码运行起来？我知道这段代码很长，可能不太容易理解，但我想做的是写一个模拟，模拟一个我之前在这里提过的问题：https://math.stackexchange.com/questions/4876146/pendulum-hanging-on-a-spinning-disk
我想做一个不错的模拟，效果应该和有人回答那个问题时提供的图片差不多。那个答案里的图片是用mathematica写的，我完全不知道怎么把它转成我需要的格式。希望你们能帮我完成这个。
代码里有两个部分。一个是计算二阶常微分方程（ODE），另一个是把结果画出来三次。当你画出这个ODE时，你会发现图线并没有按照预期的样子走。我不知道错误出在哪里，但希望你们能帮忙找出问题。
下面是这两段代码：
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import display
from sympy.vector import CoordSys3D
from scipy.integrate import solve_ivp
def FindDGL():
N = CoordSys3D('N')
e = N.i + N.j + N.k
t = sp.symbols('t')
x = sp.symbols('x')
y = sp.symbols('y')
z = sp.symbols('z')
x = sp.Function('x')(t)
y = sp.Function('y')(t)
z = sp.Function('z')(t)
p = x*N.i + y*N.j + z*N.k
m = sp.symbols('m')
g = sp.symbols('g')
r = sp.symbols('r')
omega = sp.symbols('omega')
q0 = sp.symbols('q0')
A = sp.symbols('A')
l = sp.symbols('l')
xl = sp.symbols('xl')
yl = sp.symbols('yl')
zl = sp.symbols('zl')
dpdt = sp.diff(x,t)*N.i + sp.diff(y,t)*N.j + sp.diff(z,t)*N.k
#Zwang = ((p-q)-l)
Zwang = (p.dot(N.i)**2*N.i +p.dot(N.j)**2*N.j +p.dot(N.k)**2*N.k - 2*r*(p.dot(N.i)*N.i*sp.cos(omega*t)+p.dot(N.j)*N.j*sp.sin(omega*t))-2*q0*(p.dot(N.k)*N.k) + r**2*(N.i*sp.cos(omega*t)**2+N.j*sp.sin(omega*t)**2)+q0**2*N.k) - l**2*N.i - l**2*N.j -l**2*N.k
display(Zwang)
dpdtsq = dpdt.dot(N.i)**2*N.i + dpdt.dot(N.j)**2*N.j + dpdt.dot(N.k)**2*N.k
#La = 0.5 * m * dpdtsq - m * g * (p.dot(N.k)*N.k) + (ZwangA*A)
L = 0.5 * m * dpdtsq + m * g * (p.dot(N.k)*N.k) - Zwang*A
#display(La)
display(L)
Lx = L.dot(N.i)
Ly = L.dot(N.j)
Lz = L.dot(N.k)
Elx = sp.diff(sp.diff(Lx,sp.Derivative(x,t)), t) + sp.diff(Lx,x)
Ely = sp.diff(sp.diff(Ly,sp.Derivative(y,t)), t) + sp.diff(Ly,y)
Elz = sp.diff(sp.diff(Lz,sp.Derivative(z,t)), t) + sp.diff(Lz,z)
display(Elx)
display(Ely)
display(Elz)
ZwangAV = (sp.diff(Zwang, t, 2))/2
display(ZwangAV)
ZwangA = ZwangAV.dot(N.i)+ZwangAV.dot(N.j)+ZwangAV.dot(N.k)
display(ZwangA)
Eq1 = sp.Eq(Elx,0)
Eq2 = sp.Eq(Ely,0)
Eq3 = sp.Eq(Elz,0)
Eq4 = sp.Eq(ZwangA,0)
LGS = sp.solve((Eq1,Eq2,Eq3,Eq4),(sp.Derivative(x,t,2),sp.Derivative(y,t,2),sp.Derivative(z,t,2),A))
#display(LGS)
#display(LGS[sp.Derivative(x,t,2)].free_symbols)
#display(LGS[sp.Derivative(y,t,2)].free_symbols)
#display(LGS[sp.Derivative(z,t,2)])
XS = LGS[sp.Derivative(x,t,2)]
YS = LGS[sp.Derivative(y,t,2)]
ZS = LGS[sp.Derivative(z,t,2)]
#t_span = (0, 10)
dxdt = sp.symbols('dxdt')
dydt = sp.symbols('dydt')
dzdt = sp.symbols('dzdt')
#t_eval = np.linspace(0, 10, 100)
XSL = XS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
YSL = YS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
ZSL = ZS.subs({ sp.Derivative(y,t):dydt, sp.Derivative(z,t):dzdt, sp.Derivative(x,t):dxdt, x:xl , y:yl , z:zl})
#display(ZSL.free_symbols)
XSLS = str(XSL)
YSLS = str(YSL)
ZSLS = str(ZSL)
replace = {"xl":"x","yl":"y","zl":"z","cos":"np.cos", "sin":"np.sin",}
for old, new in replace.items():
XSLS = XSLS.replace(old, new)
for old, new in replace.items():
YSLS = YSLS.replace(old, new)
for old, new in replace.items():
ZSLS = ZSLS.replace(old, new)
return[XSLS,YSLS,ZSLS]
Result = FindDGL()
print(Result[0])
print(Result[1])
print(Result[2])
这是第二段：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
def Q(t):
omega = 1
return r * (np.cos(omega * t) * np.array([1, 0, 0]) + np.sin(omega * t) * np.array([0, 1, 0])) + np.array([0, 0, q0])
def equations_of_motion(t, state, r, omega, q0, l):
x, y, z, xp, yp, zp = state
dxdt = xp
dydt = yp
dzdt = zp
dxpdt = dxdt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*r**2*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*x*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*r**2*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*q0*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*q0*x/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*x*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*x*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*y*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x**2*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
dypdt = dxdt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*r**2*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*r**2*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*y*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*q0*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*q0*y/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + g*y*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*x*np.sin(omega*t)*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*r**2*y*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*y*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*y**2*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
dzpdt = dxdt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dxdt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dxdt*omega*q0*r*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dxdt*omega*r*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dydt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dydt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - 2.0*dydt*omega*q0*r*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*dydt*omega*r*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + dzdt**2*q0/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - dzdt**2*z/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r**2*np.sin(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*r**2*np.cos(omega*t)**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*g*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + 2.0*g*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*x**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - g*y**2/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*q0*r*x*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) + omega**2*q0*r*y*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*x*z*np.cos(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2) - omega**2*r*y*z*np.sin(omega*t)/(q0**2 - 2.0*q0*z + r**2*np.sin(omega*t)**2 + r**2*np.cos(omega*t)**2 - 2.0*r*x*np.cos(omega*t) - 2.0*r*y*np.sin(omega*t) + x**2 + y**2 + z**2)
return [dxdt, dydt, dzdt, dxpdt, dypdt, dzpdt]
r = 0.5
omega = 1.2
q0 = 1
l = 1
g = 9.81
#{x[0] == r, y[0] == x'[0] == y'[0] == z'[0] == 0, z[0] == q0 - l}
initial_conditions = [r, 0, 0, 0, 0, q0-l]
tmax = 200
solp = solve_ivp(equations_of_motion, [0, tmax], initial_conditions, args=(r, omega, q0, l), dense_output=True)#, method='DOP853')
t_values = np.linspace(0, tmax, 1000)
p_values = solp.sol(t_values)
print(p_values.size)
d =0.5
Qx = [Q(ti)[0] for ti in t_values]
Qy = [Q(ti)[1] for ti in t_values]
Qz = [Q(ti)[2] for ti in t_values]
fig = plt.figure(figsize=(20, 16))
ax = fig.add_subplot(111, projection='3d')
ax.plot(p_values[0], p_values[1], p_values[2], color='blue')
ax.scatter(r, 0, q0-l, color='red')
ax.plot([0, 0], [0, 0], [0, q0], color='green')
ax.plot(Qx, Qy, Qz, color='purple')
#ax.set_xlim(-d, d)
#ax.set_ylim(-d, d)
#ax.set_zlim(-d, d)
ax.view_init(30, 45)
plt.show()
1 个回答
0
编辑：为了包含在笛卡尔坐标系中的直接解决方案，因为这正是提问者最初的方向。请查看这个回答的底部。
@Mo711，你的代码（特别是第二个示例）太长了，难以检查。拉格朗日力学的应用其实很简单——只要你选择合适的广义坐标。
这是一个有两个自由度的问题（所以你在Stack Exchange上的第一篇帖子肯定是错的：它只有一个自由度）。在我看来，使用三个笛卡尔坐标并且还要施加约束是很不自然的，因为它们并不是独立的。相反，正如Stack Exchange上某位回复者所建议的，我会使用两个角度。下面的分析会得出类似于SE回复中的角度方程，但并不完全相同：我认为他/她是从水平面取角度θ，而我下面使用的是从竖直向下的角度。
首先是拉格朗日方程。
设ϕ为包含绳子的竖直平面与x轴之间的角度。设θ为绳子与竖直向下之间的角度。那么，摆锤的坐标（相对于环的中心）为：
对时间求导，我们得到速度分量
通过求和、平方和使用三角公式简化：
拉格朗日量（严格来说，是拉格朗日量除以质量，但在后续分析中质量会抵消）为
因此
从保守系统的拉格朗日方程中
我们经过大量代数运算得到我们自由度的关键方程：
下面的代码示例使用solve_ivp来解决这些方程。请注意，第二个方程的分母意味着你不应该从摆锤竖直的位置开始，因为sinθ会是0。
然后是你的动画。
我使用FuncAnimation来自matplotlib.animation。这里的目标是更新你图表中的任何元素（在这个例子中是摆锤绳子的末端）在每一帧中。代码默认绘制动画，但你可以去掉后面的注释，将其存储为一个（相当大的）图形文件。
如果你想要轨迹而不是动画，那么在底部重新注释命令，选择plot_figure。
请注意，这个系统是混沌的，对盘子和摆锤的相对大小以及角速度ω非常敏感。如果我进行维度分析，轨迹形状应该是g/ω²L（或者盘子和独立摆锤的周期比）、R/L和初始角度的函数。
代码
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
g = 9.81
def plot_animation( t, qx, qy, qz, x, y, z ):
plotInterval = 1
fig = plt.figure( figsize=(4,4) )
ax = fig.add_subplot(111, projection='3d' )
ax.view_init(30, 45)
ax.set_xlim( -1.0, 1.0 ); ax.set_ylim( -1.0, 1.0 ); ax.set_aspect('equal')
ax.plot( qx, qy, qz, 'k' ) # ring
a = ax.plot( [qx[0],x[0]], [qy[0],y[0]], [qz[0],z[0]], 'g' ) # pendulum string
b = ax.plot( [x[0]], [y[0]], [z[0]], 'ro' ) # pendulum bob
def animate( i ):
a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] ) # update anything that has changed
b[0].set_data_3d( [x[i]], [y[i]], [z[i]] )
ani = animation.FuncAnimation( fig, animate, interval=4, frames=len( t ) )
plt.show()
# ani.save( "demo.gif", fps=50 )
def plot_figure( qx, qy, qz, x, y, z ):
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
ax.plot( qx, qy, qz, 'k' ) # disk
ax.plot( x, y, z, 'b' ) # bob trajectory
ax.plot( ( qx[-1], x[-1] ), ( qy[-1], y[-1] ), ( qz[-1], z[-1] ), 'g' ) # final line
ax.plot( x[-1], y[-1], z[-1], 'ro' ) # final bob
ax.view_init(30, 45)
ax.set_xlim( -1.0, 1.0 ); ax.set_ylim( -1.0, 1.0 ); ax.set_aspect('equal')
plt.show()
def deriv( t, Y, R, omega, L ):
theta, phi, thetaprime, phiprime = Y
ct, st = math.cos( theta ), math.sin( theta )
phase = omega * t - phi
thetaprime2 = ( phiprime ** 2 * L * st * ct + omega ** 2 * R * ct * math.cos( phase ) - g * st ) / L
phiprime2 = ( -2 * thetaprime * phiprime * L * ct + omega ** 2 * R * math.sin( phase ) ) / ( L * st )
return [ thetaprime, phiprime, thetaprime2, phiprime2 ]
R = 0.5
omega = 2.0
L = 1.0
Y0 = [ 0.01, 0.01, 0.0, 0.0 ]
period = 2 * np.pi / omega
tmax = 5 * period
solution = solve_ivp( deriv, [0, tmax], Y0, args=(R,omega,L), rtol=1.0e-6, dense_output=True )
t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )
theta = Y[0,:]
phi = Y[1,:]
# Position on disk
qx = R * np.cos( omega * t )
qy = R * np.sin( omega * t )
qz = np.zeros_like( qx )
# Trajectory
x = qx + L * np.sin( theta ) * np.cos( phi )
y = qy + L * np.sin( theta ) * np.sin( phi )
z = - L * np.cos( theta )
#plot_figure( qx, qy, qz, x, y, z )
plot_animation( t, qx, qy, qz, x, y, z )
轨迹：
如果你绝对必须在笛卡尔坐标系中解决这个问题（我强烈建议你不要这样做），那么你可以大大简化在Stack Exchange中给出的解决方案。拉格朗日量为
其中x是摆锤的位置，而"q"=(R cosωt,R sinωt,0)。将三个拉格朗日方程重写为向量形式得到
λ必须从约束方程中找到
对其进行两次求导，并注意到
结合上面的加速度表达式，得到下面的2λ表达式。因此，你的方程组为
这比Stack Exchange中给出的表达式简单得多，但可能是等效的。下面的代码展示了这一点。与第一段代码的唯一真正区别是例程deriv()，因为你现在有不同的依赖变量。
使用约束拉格朗日的另一个问题是初始条件（位置和速度）必须同时满足约束方程及其导数。也就是说，绳长是L，并且最初没有拉长。
使用笛卡尔依赖变量的代码。
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
g = 9.81
def plot_animation( t, qx, qy, qz, x, y, z ):
plotInterval = 1
fig = plt.figure( figsize=(4,4) )
ax = fig.add_subplot(111, projection='3d' )
ax.view_init(30, 45)
ax.set_xlim( -1.0, 1.0 ); ax.set_ylim( -1.0, 1.0 ); ax.set_aspect('equal')
ax.plot( qx, qy, qz, 'k' ) # ring
a = ax.plot( [qx[0],x[0]], [qy[0],y[0]], [qz[0],z[0]], 'g' ) # pendulum string
b = ax.plot( [x[0]], [y[0]], [z[0]], 'ro' ) # pendulum bob
def animate( i ):
a[0].set_data_3d( [qx[i],x[i]], [qy[i],y[i]], [qz[i],z[i]] ) # update anything that has changed
b[0].set_data_3d( [x[i]], [y[i]], [z[i]] )
ani = animation.FuncAnimation( fig, animate, interval=4, frames=len( t ) )
plt.show()
# ani.save( "demo.gif", fps=50 )
def plot_figure( qx, qy, qz, x, y, z ):
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
ax.plot( qx, qy, qz, 'k' ) # disk
ax.plot( x, y, z, 'b' ) # bob trajectory
ax.plot( ( qx[-1], x[-1] ), ( qy[-1], y[-1] ), ( qz[-1], z[-1] ), 'g' ) # final line
ax.plot( x[-1], y[-1], z[-1], 'ro' ) # final bob
ax.view_init(30, 45)
ax.set_xlim( -1.0, 1.0 ); ax.set_ylim( -1.0, 1.0 ); ax.set_aspect('equal')
plt.show()
def deriv( t, Y, R, omega, L ):
x, y, z, xdot, ydot, zdot = Y
qx, qy = R * math.cos( omega * t ), R * math.sin( omega * t )
qxdot, qydot = -omega * R * math.sin( omega * t ), omega * R * math.cos( omega * t )
twoLambda = ( g * z + ( R ** 2 - qx * x - qy * y ) * omega ** 2 - ( xdot - qxdot ) ** 2 - ( ydot -qydot ) ** 2 - zdot ** 2 ) / L ** 2
xddot = twoLambda * ( x - qx )
yddot = twoLambda * ( y - qy )
zddot = twoLambda * z - g
return [ xdot, ydot, zdot, xddot, yddot, zddot ]
R = 0.5
omega = 2.0
L = 1.0
Y0 = [ R, 0.0, -L, 0.0, 0.0, 0.0 ]
period = 2 * np.pi / omega
tmax = 5 * period
solution = solve_ivp( deriv, [0, tmax], Y0, args=(R,omega,L), rtol=1.0e-6, dense_output=True )
t = np.linspace( 0, tmax, 1000 )
Y = solution.sol( t )
x = Y[0,:]
y = Y[1,:]
z = Y[2,:]
# Position on disk
qx = R * np.cos( omega * t )
qy = R * np.sin( omega * t )
qz = np.zeros_like( qx )
#plot_figure( qx, qy, qz, x, y, z )
plot_animation( t, qx, qy, qz, x, y, z )