计算科学 - 已解决：阻尼弹簧质量系统、错误的位置、正确的速度和加速度 - 吾爱随笔录

我正在使用 python 编程调制具有重力和空气阻力的弹簧质量系统。弹簧垂直悬挂并附有重物。然后用户选择一个长度将其向下拖动，代码将生成图形

我用这种方法解了微分方程。

\begin{aligned} \sum F_{x} & = 0 \\ \sum F_{y} & = - k y + m g - F_{drag} \\ F & = m a \\ m a_{y} & = - k y + m g + F_{drag} \\ m \frac{d v}{d t} & = - k y + m g - F_{drag} \\ m d v & = - k y d t + m g d t - F_{drag} d t \\ m Δ v & = - k y Δ t + m g Δ t - F_{drag} Δ t \\ Δ t & = 0.001 s \\ Δ v & = \frac{\sum F_{y}}{m} \\ v^{i + 1} & = v^{i} + Δ v \\ Δ s & = v Δ t \\ x & = x + Δ s \end{aligned}

$\begin{align} \sum F_x &= 0\\ \sum F_y &= -ky + mg - F_\text{drag}\\ \mathbf{F} &= m \mathbf{a}\\ m a_y &= -ky + mg + F_\text{drag}\\ m\frac{\mathrm{d}v}{\mathrm{d}t} &= -ky + mg - F_\text{drag}\\ m \mathrm{d}v &= - ky \mathrm{d}t + mg \mathrm{d}t - F_\text{drag}\mathrm{d}t\\ m\Delta v &= -ky\Delta t + mg\Delta t - F_\text{drag}\Delta t\\ \Delta t &= 0.001 \text{ s}\\ \Delta v &= \frac{\sum F_y}{m}\\ v^{i + 1} &= v^{i} + \Delta v\\ \Delta s &= v \Delta t\\ x &= x + \Delta s \end{align}$

我用这段代码在python中解决了这个问题：

from pylab import *
from matplotlib import pyplot as plt # kilde: https://matplotlib.org/users/pyplot_tutorial.html

# housekeeping
g = 9.81  # gravitation
m = float(input("how much do you want the mass to be? "))
l = 0.02     # aerodrag constant
k = float(input("insert springconstant:"))
#startvalues
t0 = 0.0  # start time
v0 = 0.0   # start speed
dv0 = 0.0#start axcelleration
xstart = float(input("how long do you want do drag it down? "))
x0 = (-m*g/k) - xstart  # start position 



time = list() # list for time
time.append(t0) # insert start value for time
speed= list() 
speed.append(v0) 
position = list() 
position.append(x0 + mg/k) 
aksellerasjon = list() # accelleration
aksellerasjon.append(dv0) 
dt= 0.001 
x = x0  
v = v0
t = t0


# main algorithm 
while t < 20.0:  

    aerodrag= l*v
    F = (-m*g - k*(x) - aerodrag)  #-m*g because gravitation is negative    
    dv = F*dt/m
    v+= dv 
    dx=  v*dt
    x += (dx) 
    t += dt
    aksellerasjon.append(dv) 
    time.append(t)
    speed.append(v)
    position.append(x)

# Plotting  
plt.figure(figsize=(6,10)) 
subplot(3,1,1)                  
title("Akselerasjon") 
xlabel("tid/s")   
ylabel("akselerasjon m/s^2")  
plot(tid, aksellerasjon, "b-")

subplot(3, 1, 2)
plot(tid, fart, "r-")
title("Hastighet")  
xlabel("tid/s")   
ylabel("fart m/s")    

subplot(3,1,3)
plot(tid, posisjon, "g-")
title("Strekning")              
xlabel("tid/s")              
ylabel("posisjon/ m")      
plt.tight_layout() 

show()

我已经用我的语言翻译了注释和部分代码，很抱歉我的英语不好

速度和加速度是正确的，但位置是错误的。它太低了，尽管波接缝的幅度正确。

我在 GeoGebra 中求解了微分方程，位置应该是这样的。

我究竟做错了什么？

在评论中得到很好的帮助后，我在附录中添加了 + mg/k

posisjon.append(x + (m*g)/k)

现在位置收敛到零，谢谢！