0
$\begingroup$

I am stepping in the coding of motion simulation and have found in internet following python codes for practice purpose.

# Simulation of Mass-Spring-Damper System
import numpy as np
import matplotlib.pyplot as plt
# Model Parameters
c = 4 # Damping constant
k = 2 # Stiffness of the spring
m = 20 # Mass
F = 5 # Force
# Simulation Parameters
Ts = 0.1
Tstart = 0
Tstop = 60
N = int((Tstop-Tstart)/Ts) # Simulation length
x1 = np.zeros(N+2)
x2 = np.zeros(N+2)
x1[0] = 0 # Initial Position
x2[0] = 0 # Initial Speed
a11 = 1
a12 = Ts
a21 = -(Ts*k)/m
a22 = 1 - (Ts*c)/m
b1 = 0
b2 = Ts/m
# Simulation
for k in range(N+1):
    x1[k+1] = a11 * x1[k] + a12 * x2[k] + b1 * F
    x2[k+1] = a21 * x1[k] + a22 * x2[k] + b2 * F

# Plot the Simulation Results
t = np.arange(Tstart,Tstop+2*Ts,Ts)
#plt.plot(t, x1, t, x2)
plt.plot(t,x1)
plt.plot(t,x2)
plt.title('Simulation of Mass-Spring-Damper System')
plt.xlabel('t [s]')
plt.ylabel('x(t)')
plt.grid()
plt.legend(["x1", "x2"])
plt.show()

with the output plot

enter image description here

  1. By setting the damping coefficient C=0, I expected the amplitudes of the object's oscillation curve stays unchanged. But the curve shows increased amplitude over time. Does anyone know why the undamped object (C=0) has the oscillation with increased amplitude?

enter image description here

  1. The increment / time step seems to have influence on oscillation curve: by changing from 0.1 to 0.001 the curve shows slight change, while by increasing from 0.1 to 2 the curve becomes diverging (numerical instability - my assumption) Can anyone share your opinion/experiences herewith? Any rules/instructions for selecting the time step for reliable results?

enter image description here

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ First realize that you are doing numerical integration with the most rudimentary method of assuming your integrant is a rectangular. What your seeing is an integration error. So first order of business is to choose a better approximation method i suggest 4th order runge-kutta its easy to implement, but you could use even just trapezoidal method (if your purpose is learning), if you write it to a function you can just reuse the function with no change to code. On the otherhand numpy has orders of magnitudes better solvers available too. Then you need to consider doing implicit solver $\endgroup$
    – joojaa
    Commented Aug 29 at 7:23
  • 1
    $\begingroup$ After this you worry about step size maybe make a adaptive solver. But really use a ready made solver if your not doing this for pure learning. $\endgroup$
    – joojaa
    Commented Aug 29 at 7:24
  • 1
    $\begingroup$ I think what @joojaa suggested is basically Crank-Nicholson method. I have used it in the past for similar problem. $\endgroup$
    – Tomáš Létal
    Commented Aug 29 at 7:51
  • 1
    $\begingroup$ Many thanks for the information and suggestion. @joojaa : after doing reading in internet about the suggested terms, i understand the your point about the the rudimentary method. I understand the suggested runge-kutta 4th order is a fixed-step and classical method, and the suggested adaptive solver (e.g. Scipy Odeint) is more advanced and modern method. Is my understanding correct? $\endgroup$
    – Daniel Liu
    Commented Aug 29 at 11:19
  • $\begingroup$ A very related video. $\endgroup$
    – AJN
    Commented Aug 29 at 12:04

1 Answer 1

Reset to default
0
$\begingroup$

It looks to me like a larger time step leads to lower accuracy. If you set the time step very small, and damping to 0, you do get what looks like a constant amplitude. With c = 0, Ts = 0.00005 (and Tstop = 240, to get more time for divergence to build up, if there is any) I get this plot:

enter image description here

So yes, c=0 does produce the un-damped oscillation you expect, if you simulate with fine enough time steps.

Oh, I should probably add that with such a small time step, and large range, this does run pretty slowly.

Improve this answer
$\endgroup$
1
  • $\begingroup$ many thanks for the explanation! $\endgroup$
    – Daniel Liu
    Commented Aug 29 at 10:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.