1

I'm trying to plot the function m*x''(t) + r*x'(t)+k*x(t) = 0 and so after some solving it goes to:

x(t+2h) = x(t+h)*(2-(r*h/m)) + x(t) * (-1 + (r*h/m) - (k*h*h/m))

so i tried making a python script that makes this recursion and then plots the graph but i'm doing sth wrong

import matplotlib.pyplot as plt
m = 0.015
k = 5.0
r = 0.24
h = 0.1
x0 = 0.05
v0 = 0


def displacement(t):
    if t <= 0:
        return x0
    elif t == h:
        return v0*h + x0
    else:
        return displacement(t-h)*(2-((r*h)/m)) + displacement(t-2*h)*(((r*h)/m) - 1 - ((k*h*h)/m))

#xn(2 - s*h/m - k*h^2/m) + xn-1(s*h/m - 1) = xn+1
# x(n - 1) = x(t - h)
#x1' = x1-x0/h => x1 = v0*h + x0

t_values = []
for i in range(50):
    t_values.append(i/10)

print(t_values)

xt_values = []

for j in t_values:
    xt_values.append(displacement(j))
    print(xt_values)

print(xt_values)

#t_values = [t/100000 for t in range(200000)]
#xt_values = [displacement(t) for t in t_values]

plt
plt.plot(t_values, xt_values)
plt.title('Damped Spring System')
plt.xlabel('Time')
plt.ylabel('Displacement(x)')
plt.show()

All things m, k, r, h are parameters as well as the first 2 cases of recursion which are 1 and 0. Is there any coding mistake I'm making or is it just a skill issue in the physics part?

The expected graph is something like a COS graph but the amplitude decreases over time and the period gets larger and larger. The graph I'm getting is a very strange alternating + and - graph that gets bigger and bigger over time almost exponentially (very fast, the computer can't handle samples very quickly.)

Smasher
  • 11
  • 2
  • 1
    You haven't said what **the actual problem** is. What does the script actually do, and what did you expect instead? – John Gordon Apr 30 '23 at 20:02
  • @JhonGordon You are right, I have edited it to the best of my abilities – Smasher Apr 30 '23 at 20:09
  • not sure this is necessarily the cause of your issue, but `elif t == h` may be a problem, because of floating point rounding issues. It might be better to do `elif abs(t - h) < some_suitably_small_epsilon:` – slothrop Apr 30 '23 at 20:10
  • It probably won't become apparent for a simulation with only 50 time steps, but pure recursion is a very inefficient approach, because every previous value of `displacement` will be recalculated on every time step. When you have the code working for a small number of steps, look at **memoization** to enable it to work for a longer simulation: https://stackoverflow.com/questions/1988804/what-is-memoization-and-how-can-i-use-it-in-python – slothrop Apr 30 '23 at 20:26
  • @slothrop Yep i tried it with no recursion but just setting 3 variables and adding to a list to it be it slower but just to check my work. It's not even close to what I want I to be but physics side also look good so I'm a bit confused – Smasher May 01 '23 at 00:00
  • Does it work any better with a smaller timestep? Generally a method like this relies on the timestep being very small (in the context of the system). Here, `ω = √(k/m) ~ 18.26`, and the period of the oscillator (ignoring damping for now) is `2π/ω ~ 0.344`. Your timestep `h = 0.1` is rather large relative to the system - the system performs about a third of a full cycle of oscillation during every timestep - so, loosely speaking, you might be "missing" behaviour by chunking time into overly large intervals. – slothrop May 01 '23 at 09:11

0 Answers0