Welcome, guest | Sign In | My Account | Store | Cart

Spring-Mass System Simulation (Python recipe) by FB36
ActiveState Code (http://code.activestate.com/recipes/577681/)

It simulates a damped spring-mass system driven by sinusoidal force.

Python, 81 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
# Damped spring-mass system driven by sinusoidal force
# FB - 201105017
import math
from PIL import Image, ImageDraw
imgx = 800
imgy = 600
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)

# Second Order ODE (y'' = f(x, y, y')) Solver using Euler method
# n : number of steps (higher the better)
# xa: initial value of independent variable
# xb: final value of independent variable
# ya: initial value of dependent variable
# y1a: initial value of first derivative of dependent variable
# Returns value of y, y1 at xb. 
def Euler2(f, xa, xb, ya, y1a, n):
      h = (xb - xa) / float(n)
      x = xa
      y = ya
      y1 = y1a
      for i in range(n):
          y1 += h * f(x, y, y1)
          y += h * y1
          x += h
      return [y, y1]

# Damped spring-mass system driven by sinusoidal force
# y'' = (F0 * math.cos(omega * t - phi) - b * y' - k * y) / m
# y'' : acceleration
# y' : velocity
# y : position
m = 2.0 # mass (kg)
F0 = 4.76 # force amplitude constant (N)
omega = 0.36 # angular frequency (rad/s)
phi = 0.0 # phase constant (rad)
b = 0.0 # friction constant (Ns/m)
k = 20.0 # spring constant (N/m)
def f(x, y, y1):
    return (F0 * math.cos(omega * x - phi) - b * y1 - k * y) / m

yaSim = 0.0 # initial position (m)
y1aSim = 0.0 # initial velocity (m/s)
n = 1000 # number of steps for Euler method
xaSim = 0.0 # initial time of simulation (s)
xbSim = 100.0 # final time of simulation (s)
xdSim = xbSim - xaSim # deltaT of simulation
nSim = 1000 # number of time steps of simulation

# find min and max values of position (needed for the graph)
ya = yaSim
y1a = y1aSim
yMin = ya
yMax = ya
for i in range(nSim):
    xa = xaSim + xdSim * i / nSim
    xb = xaSim + xdSim * (i + 1) / nSim
    y_y1 = Euler2(f, xa, xb, ya, y1a, n)
    ya = y_y1[0]
    y1a = y_y1[1]
    if ya < yMin:
        yMin = ya
    if ya > yMax:
        yMax = ya

# draw the graph
ya = yaSim
y1a = y1aSim
for i in range(nSim):
    xa = xaSim + xdSim * i / nSim
    xb = xaSim + xdSim * (i + 1) / nSim
    kxa = (imgx - 1) * (xa - xaSim) / xdSim
    kya = (imgy - 1) * (ya - yMin) / (yMax - yMin)
    y_y1 = Euler2(f, xa, xb, ya, y1a, n)
    ya = y_y1[0]
    y1a = y_y1[1]
    kxb = (imgx - 1) * (xb - xaSim) / xdSim
    kyb = (imgy - 1) * (ya - yMin) / (yMax - yMin)
    draw.line((kxa, kya, kxb, kyb), (0, 255, 0)) # (r, g, b)

image.save("Spring_mass system simulation.png", "PNG")
Tags: chaos, math, mathematics, physics, simulation

1 comment

FB36 (author) 12 years, 11 months ago  # | flag

The spring-mass system equation section can be replaced to simulate a driven pendulum:

  # Forced Damped Pendulum
  # f(t)-b*l*y'-m*g*sin(y)=m*l*y'' and f(t)=F0*cos(omega*t-phi) then
  # y''=(F0*cos(omega*t-phi)-b*l*y'-m*g*sin(y))/(m*l)
  # y'' : acceleration
  # y' : velocity
  # y : position
  m = 1.0 # mass (kg)
  F0 = 4.76 # force amplitude constant (N)
  omega = 0.36 # angular frequency (rad/s)
  phi = 0.0 # phase constant (rad)
  b = 0.0 # friction constant (Ns/m)
  g = 9.807 # standard gravity (m/s^2)
  L = 1.0 # length of pendulum (m)
  def f(x, y, y1):
      return (F0*math.cos(omega*x-phi)-b*L*y1-m*g*math.sin(y))/(m*L)