弹性摆系统的runge-kutta

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle l0 = 10 #spring at rest g = 9.81 #gravity m = 1 #mass of particle k = 40 #spring constant dt = 0.1 #upgrade Theta0 = 3*np.pi/4 #initial theta z10 = 0 #initial theta velocity z20 = 0 #initial l velocity tmax, dt = 20, 0.01 t = np.arange(0, tmax+dt, dt) def f_theta(z1, z2, theta, g, L): return (-g*np.sin(theta) - 2*z1*z2) / L def f_L(z1,theta, g, L, l0, m, k): return (m*L*z1**2 - k*(L-l0) + m*g*np.cos(theta)) / m Thetapoints = [] z1 = [] Lpoints = [] z2 = [] for x in t: Thetapoints.append(Theta0) z1.append(z10) Lpoints.append(l0) z2.append(z20) m1 = dt*z10 M1 = dt*f_theta(z10,z20,Theta0,g,l0) k1 = dt*z20 K1 = dt*f_L(z10,Theta0,g,l0,l0,m,k) m2 = dt*(z10+0.5*M1) M2 = dt*(f_theta(z10+0.5*M1,z20+0.5*K1,Theta0+0.5*m1,g,l0+0.5*k1)) k2 = dt*(z20+0.5*K1) K2 = dt*(f_L(z10+0.5*M2,Theta0+0.5*m2,g,l0+0.5*k2,l0,m,k)) m3 = dt*(z10+0.5*M2) M3 = dt*f_theta(z10+0.5*M2,z20+0.5*K2,Theta0+0.5*m2,g,l0+0.5*k2) k3 = dt*(z20+0.5*K2) K3 = dt*(f_L(z10+0.5*M2,Theta0+0.5*m2,g,l0+0.5*k2,l0,m,k)) m4 = dt*(z10+M3) M4 = dt*f_theta(z10+M3,z20+K3,Theta0+m3,g,l0+k3) k4 = dt*(z20+K3) K4 = dt*(f_L(z10+M3,Theta0+m3,g,l0+k3,l0,m,k)) Theta0 += (m1 + 2*m2 + 2*m3 + m4)/6 l0 += (k1 + 2*k2 + 2*k3 + k4)/6 z10 += (M1 + 2*M2 + 2*M3 + M4)/6 z20 += (K1 + 2*K2 + 2*K3 + K4)/6 x = np.array(Lpoints)* np.sin(np.array(Thetapoints)) y = -np.array(Lpoints)* np.cos(np.array(Thetapoints)) plt.plot(t,Lpoints) plt.show()

1条回答

网友

1楼 · 发布于 2024-09-27 07:23:56

最明显的错误是l0不仅用作弹簧的恒定静止长度，还用作弹簧的动态长度，结果不可预测

在更系统的方法中，将系统编码为系统，并使用矢量版本的RK4

def federpendel(u,m,g,l0,k):
    th, r, Vth, Vr = u
    Ath = (-g*np.sin(th) - 2*Vth*Vr) / r
    Ar  = r*Vth**2 - k/m*(r-l0) + g*np.cos(th)
    return np.array([ Vth, Vr, Ath, Ar])

l0 = 10             #spring at rest
g  = 9.81          #gravity
m  = 1             #mass of particle
k  = 40            #spring constant
th0 = 3*np.pi/4 #initial theta
Dth0    = 0         #initial theta velocity
Dr0    = 0         #initial  l velocity
u = np.array([ th0, l0, Dth0, Dr0])
dt = 0.1           #upgrade 
tmax= 20
t = np.arange(0, tmax+0.5*dt, dt)
U = [u];
for n in range(len(t)-1):
    k1 = federpendel(u,m,g,l0,k)*dt
    k2 = federpendel(u+0.5*k1,m,g,l0,k)*dt
    k3 = federpendel(u+0.5*k2,m,g,l0,k)*dt
    k4 = federpendel(u+k3,m,g,l0,k)*dt
    u = u + (k1+2*k2+2*k3+k4)/6
    U.append(u)

th, r, Dth, Dr = np.asarray(U).T
plt.subplot(2,1,1);plt.plot(t,r);
plt.subplot(2,1,2);plt.plot(t,th);
plt.show()

给出情节

相关问题更多 >

编程相关推荐

热门问题

热门文章