Моделирование маятника, висящего на вращающемся диске
Может ли кто-нибудь запустить этот код? Я знаю, что это очень долго и, возможно, нелегко понять, но я пытаюсь написать симуляцию проблемы, которую я уже публиковал здесь: https://math.stackexchange.com/ вопросы/4876146/маятник-висит на вращающемся диске
Я пытаюсь создать красивую симуляцию, которая бы выглядела так, как будто кто-то ответил на связанный вопрос. Картинка в ответе написана на языке математики, и я понятия не имел, как ее перевести. Надеюсь, ты поможешь мне закончить это.
Есть два элемента кода. Один рассчитывает ОДУ второй степени, а другой строит его 3 раза. Когда вы строите ОДУ, вы можете видеть, что линия графика не делает того, что должна. Я не знаю, где ошибка, но надеюсь, вы сможете помочь.
Вот два фрагмента:
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()
Извините, что предоставил всю свою кодовую базу, но вопрос наполовину о ODE, который он создает, наполовину о построении графиков для визуализации. Где недостающее звено между правильным лагранжианом и неправильным моделированием в трехмерном пространстве. Код работает, и вопрос в том, почему симуляция выглядит неправильно. Возможно, этот вопрос относится к обмену стеками математики, но я хотел попробовать спросить людей здесь.
В нынешнем виде @ Mo711 ваш код настолько сложен, что его невозможно проверить. Неясно даже, верны ли лежащие в основе уравнения. В ссылке на Stack Exchange один ответ включает гораздо лучшее предложение использовать два полярных угла. Лучше использовать 2 степени свободы, чем 3 обобщенные координаты с отдельным ограничением. (Обратите внимание, что вы можете захотеть переработать решение с двумя углами, поскольку используемые углы на самом деле не те, которые предложены.) При таком подходе с двумя углами его легко закодировать на Python. Но, к сожалению, тема закрыта.
@lastchance, давайте найдём ещё одного человека, который проголосует и возобновит этот вопрос!
@Davide_sd, да, было бы хорошо снова открыть его: это была увлекательная проблема, и анимация на Python хороша. К сожалению, мне все еще не хватает репутации, чтобы проголосовать за возобновление работы.
@Mike'Pomax'Kamermans, пожалуйста, нужен еще один голос, чтобы снова открыть этот вопрос. Этот вопрос важен, потому что модуль механики Sympy очень мощный, но очень сложен в использовании. На самом деле, многие вопросы по этой теме так и не получают ответа. На данный момент потенциально есть два человека, которые могут предложить решение этого вопроса, что определенно даст много информации о симуляциях с Sympy.
очень хорошо, теперь вы сможете написать ответ. Но этот пост все еще нуждается в улучшении потому что, как написано, вопрос в том, «может ли кто-нибудь заставить это запустить», а это не то, для чего нужен SO. @ Mo711, пожалуйста, немного обновите сообщение, чтобы объяснить реальную проблему с кодом, с которой вы столкнулись, в каких строках, по вашему мнению, что-то идет не так, и где, основываясь на том, какую отладку и (повторный) поиск вы уже выполнили, потому что вы пишу не для себя. Вы публикуете сообщение для себя и других людей, у которых возникнет такая же проблема в будущем, и они должны иметь возможность сказать, что у вас возникла их проблема.
Ответы 1
Отредактировано, чтобы включить прямое решение в декартианах, поскольку это было первоначальное направление, выбранное ОП. См. нижнюю часть этого ответа.
Снова отредактировано, чтобы обеспечить альтернативный вывод декартовых уравнений из Второго закона Ньютона (F = ma), а не из лагранжевой механики. (Дополнительным преимуществом этого метода является то, что он также определяет натяжение стержня маятника.)
Уравнения Лагранжа (i) Угловые координаты
Этот метод, связанный с ответами в вашем сообщении Stack Exchange, но не совсем такой же, весьма привлекателен, поскольку он решает вашу проблему с двумя степенями свободы с двумя независимыми координатами, поэтому не приходится иметь дело с ограничением (которое длина маятника фиксирована).
Пусть φ — угол между вертикальной плоскостью, содержащей струну, и осью x. Пусть θ — угол, образуемый струной с нисходящей вертикалью. Тогда координаты боба (относительно центра кольца):
Дифференцируем по времени и получаем компоненты скорости
Суммирование, возведение в квадрат и упрощение с использованием тригонометрических формул:
Лагранжиан (строго говоря, лагранжиан, разделенный на массу, но в последующем анализе масса будет сокращаться) равен
откуда
Из уравнений Лагранжа для консервативной системы
мы получаем (после МНОГО алгебры!) ключевые уравнения для наших степеней свободы:
В приведенном ниже примере кода эти проблемы решаются с помощью solve_ivp. Обратите внимание, что знаменатель второго уравнения означает, что вам не следует начинать с вертикали маятника, поскольку sinθ будет равна 0.
Затем ваша анимация.
Я использую FuncAnimation из matplotlib.animation. Целью здесь является обновление любых элементов вашего графика (в данном случае концов струны маятника) в каждом кадре. По умолчанию код отображает анимацию, но вы можете удалить последующий комментарий и сохранить его как (довольно большой) графический файл.
Если вам нужна траектория, а не анимация, повторно прокомментируйте команды внизу, чтобы вместо этого выбрать plot_figure.
Обратите внимание, что система хаотична и очень чувствительна к относительным размерам диска и маятника, а также угловой скорости ω. Если я проведу анализ размеров, форма траектории должна быть функцией g/ω^2 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 )
Траектория:
Уравнения Лагранжа (ii) Декартовы координаты
Трудность с декартовыми координатами заключается в том, что вы используете три координаты для решения задачи с двумя степенями свободы, и поэтому ограничение (фиксированная длина маятника) должно решаться с помощью множителя Лагранжа λ.
Следующее значительно упрощает решение, представленное в 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 )
Ньютоновская механика
Ваш недавний пост (об упругих маятниках) побудил меня спросить, как проблема может быть решена с помощью механики Ньютона F = ma, а не с помощью лагранжианов.
На качающийся маятник действуют две силы: (неизвестное) натяжение T вдоль вектора от качания до точки привязи и вес качания, мг. Разделив на массу, получим ускорение.
Нам нужно найти натяжение T. Мы делаем это, применяя ограничение, согласно которому маятник имеет фиксированную длину:
Дважды продифференцировав и разделив на 2, получим
Умножая, используя приведенное выше выражение для ускорения и отмечая, что
дает
откуда мы получаем выражение для напряжения Т:
Наконец, подставив это обратно в уравнение ускорения (и перевернув знаки для умножения на x-q, а не на q-x), после небольшой перестановки получим:
что совпадает с уравнением, полученным с помощью лагранжевой механики выше.
Большое спасибо. Мне все еще интересно, почему все это не зависит от массы, но, думаю, анализ превосходит интуицию...
Вы можете оставить массу, если хотите… но в конечном итоге она просто исчезнет, потому что она одинаково проявляется в кинетической и гравитационной потенциальной энергии. Все было бы по-другому, если бы ваша струна была эластичной, потому что она не проявлялась бы в упругой потенциальной энергии.
Извините, что так поздно спросил, но верёвку заменил на пружину. Имеет ли смысл заменить L на (L+s) для отклонения, а для потенциальной энергии просто иметь cos(phi)⋅g⋅L+k⋅s, где k — постоянный коэффициент пружины?
Если «веревка» становится эластичной, вам нужна другая степень свободы (например, расширение. Потенциальная энергия упругости определяется выражением (1/2)kx^2.
можете ли вы рассказать об этом подробнее?
В случае, когда струна эластична, у вас больше нет ограничений на длину струны. В этом случае, вероятно, проще использовать механику Ньютона (F=ma), чем лагранжево движение. Я ответил на ваш другой пост по этому поводу.
Stackoverflow предназначен не для того, чтобы помочь вам запустить всю вашу кодовую базу. Если написанный вами код не работает (или работает не так, как ожидалось) у вас, то вам придется углубиться в то, где в коде что-то впервые начинает идти не так, основываясь на знании того, что код должен делать с учетом какой-то ввод, и то, что вы начинаете видеть вместо этого (сразу или через некоторое время), а затем сосредоточьте свое исследование на этой части кода, потому что в одной строке кода что-то начинает идти не так: найдите эту строку. Даже если в последующих строках тоже могут быть проблемы. Вы найдете их дальше.