Mathematical Modelling of Spring-Mass system with Scipy

Let’s take a look at one of the problems provided by the SciPy community. The example is available on the Internet, at the SciPy website. However, some of the methods explained in this example are deprecated; hence, we’ll rebuild the example, so that it works correctly with the latest version of SciPy and NumPy.

We’re going to build and simulate a model based on a coupled spring-mass system, which is essentially a harmonic oscillator, in which a spring is stretched or compressed by a mass, thereby developing a restoring force in the spring, which results in harmonic motions when the mass is displaced from its equilibrium position.

For an undamped system, the motion of Block 1 is given by the following differential equation:

m1d2x1dt=$m1\frac{d^{2}x1}{dt}=$

(k1+k)x1−k2x2=0$(k1+k)x1-k2x2=0$

For Block 2:

m2d2x2dt$m2\frac{d^{2}x2}{dt}$

+kx2−k1x1=0$+kx2-k1x1=0$

In this example, we’ve taken a coupled spring-mass system, which is subjected to a frictional force, thereby resulting in damping. Note that damping tends to reduce the amplitude of oscillations in an oscillatory system. For our example, let us assume that the lengths of the springs, when subjected to no external forces, are L1 and L2. The following differential equations define such a system:

m1d2x1dt+μ1dx1dt+k1(x1−L1)$m1\frac{d^{2}x1}{dt}+\mu 1\frac{dx1}{dt}+k1(x1-L1)$

−k2(x1−x2−L2)=0$-k2(x1-x2-L2)=0$

and

m2d2x2dt+μ2dxdt+$m2\frac{d^{2}x2}{dt}+\mu 2\frac{dx}{dt}+$

k(x2−x1−L2)=0$k(x2-x1-L2)=0$

We’ll be using the Scipy odeint function to solve this problem. The function works for first-order differential equations; hence, we’ll re-write the equations as first fourth order equations:

dx1dt=y1$\frac{dx1}{dt}=y1$ dy1dt=(−μ1y1−k1(x1−L1)+k(x2−x1−L2)m1$\frac{dy1}{dt}$

$\; \; \; \; \; \; =\frac{(-\mu 1y1-k1(x1-L1)+k(x2-x1-L2)}{m1}$ dx2dt=y$\frac{dx2}{dt}=y$

dy2dt=(−μ2y2−k2(x2−x1−L1))m2$\frac{dy2}{dt}=\frac{(-\mu 2y2-k2(x2-x1-L1))}{m2}$

Now, let’s write a simple Python script to define this problem:

def vector(w,t,p):

x1,y1,x2,y2 = w

m1,m2,k1,k2,u1,u2,L1,L2 = p

f = [y1, (-b1y1 - k1(x1-L1) + k2*(x2-x1-L2))/m1,

y2, (-b2y2 - k2(x2-x1-L2))]

return f

In this script, we have simply defined the above mentioned equations programmatically. The argument w defines the state variables; t is for time, and p defines the vector of the parameters. In short, we have simply defined the vector field for the spring-mass system in this script.

Now, let’s define a script that uses odeint to solve the equations for a given set of parameter values, initial conditions, and time intervals. The script prints the points in the solution to the terminal.

from scipy.integrate import odeint

import two_springs

Parameter values

Masses:

m1 = 1.0

m2 = 1.5

Spring constants

k1 = 8.0

k2 = 40.0

Natural lengths

L1 = 0.5

L2 = 1.0

Friction coefficients

b1 = 0.8

b2 = 0.5

Initial conditions

x1 and x2 are the initial displacements; y1 and y2 are the initial velocities

x1 = 0.5

y1 = 0.0

x2 = 2.25

y2 = 0.0

ODE solver parameters

abserr = 1.0e-8

relerr = 1.0e-6

stoptime = 10.0

numpoints = 250

Create the time samples for the output of the ODE solver.

t = [stoptime*float(i)/(numpoints-1) for i in range(numpoints)]

Pack up the parameters and initial conditions:

p = [m1,m2,k1,k2,L1,L2,b1,b2]

w0 = [x1,y1,x2,y2]

Call the ODE solver.

wsol = odeint(two_springs.vectorfield,w0,t,args=(p,),atol=abserr,rtol=relerr)

Print the solution.

for t1,w1 in zip(t,wsol):

print t1,w1[0],w1[1],w1[2],w1[3]

The scipy.integrate.odeint function integrates a system of ordinary differential equations. It takes the following parameters:

func: callable(y, t0, ...) — It computes the derivative of y at t0. y0: array — This is the initial condition on y; it can be a vector. t: array — It is a sequence of time points for which to solve for y. The initial value point should be the first element of this sequence. args: tuple — Indicates extra arguments to pass to function. In our example, we have added atrol and rtol as extra arguments to deal with absolute and relative errors. The zip function takes one or more sequences as arguments, and returns a series of tuples that pair up parallel items taken from those sequences.

Copy the solution generated from this script to a text file using the cat command. Name this text file as two_springs.txt.

The following script uses Matplotlib to plot the solution generated by two_springs_solver.py:

from pylab import *

from matplotlib.font_manager import FontProperties

import numpy as np

data = np.loadtxt('two_springs.txt')

t = data[:,0]

x1 = data[:,1]

y1 = data[:,2]

x2 = data[:,3]

y2 = data[:,4]

figure(1,figsize=(6,4))

xlabel('t')4

grid(True)

hold(True)

lw = 1

plot(t,x1,'b',linewidth=lw)

plot(t,x2,'g',linewidth=lw)

legend((r'x1$x_1$',r'x2$x_2$'),prop=FontProperties(size=16))

title('Mass Displacements for the Coupled Spring-Mass System')

savefig('two_springs.png',dpi=72)

On running the script, we get the plot shown in Figure 12. It clearly shows how the mass displacements are reduced with time for damped systems.

I hope this article gives you an idea about Mathematical Modeliing with Scipy. We'll take up some more complex problems in forthcoming articles.