1
$\begingroup$

I am modulating a spring-mass system with gravitation and aero drag, with python programming. The spring is hanging vertically and attached a weight. The user then selects a length to drag it down and the code will generate the graphs

I solved the differential equation this way.

\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}

I solved this in python with this code:

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()

I have translated the comments and parts of the code from my language, so sorry for my bad English

The speed and acceleration are correct, but the position is wrong. It's too low, although the amplitude of the wave seams correctly.

enter image description here

I solved the differential equation in GeoGebra and the position should look like this.

enter image description here

What am I doing wrong?

After getting good help in the comments, i added + mg/k in the append

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

The position now converges to zero, thanks!

enter image description here

$\endgroup$
  • 3
    $\begingroup$ You just seem to have an offset. The damped oscillation doesn't seem to converge to a steady state at $x=0$ over time... $\endgroup$ – Wolfgang Bangerth Jun 17 at 18:33
  • $\begingroup$ Yeah, i areee, how can i adjust for this offset in the plotting? I've tried adding x start in the x parameter, but it turns out wrong. $\endgroup$ – Kathiravann Jun 18 at 15:33
  • $\begingroup$ I suggest that you edit your equations properly, so we can understand better what you are trying to do. $\endgroup$ – nicoguaro Jun 18 at 15:57
  • $\begingroup$ Hi, i have attached an image from word which better explains the mathematical method $\endgroup$ – Kathiravann Jun 18 at 16:09
  • 1
    $\begingroup$ I have typed your equations. In the future, please use MathJax for your equations. $\endgroup$ – nicoguaro Jun 18 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.