耦合弹簧-质量系统

日期2018-02-17（最后修改），2009-01-21（创建）

此 Cookbook 示例展示了如何求解微分方程组。（其他示例包括 Lotka-Volterra 教程，僵尸末日 和 KdV 示例。）

耦合弹簧-质量系统

此图显示了要建模的系统

两个质量分别为 $m_1$ 和 $m_2$ 的物体通过弹簧常数分别为 $k_1$ 和 $k_2$ 的弹簧耦合。左侧弹簧的左端固定。我们假设弹簧在不受外力作用时的长度分别为 $L_1$ 和 $L_2$。

质量在产生摩擦的表面上滑动，因此存在两个摩擦系数，$b_1$ 和 $b_2$。

该系统的微分方程为

$m_1 x_1'' + b_1 x_1' + k_1 (x_1 - L_1) - k_2 (x_2 - x_1 - L_2) = 0$

$m_2 x_2'' + b_2 x_2' + k_2 (x_2 - x_1 - L_2) = 0$

这是一个耦合的二阶方程组。为了用 SciPy 提供的 ODE 求解器之一求解该系统，我们必须首先将其转换为一阶微分方程组。我们引入两个变量

$y_1 = x_1'$

$y_2 = x_2'$

这些是质量的速度。

经过一些代数运算，我们可以将两个二阶方程改写为以下四个一阶方程组

$x_1' = y_1$

$y_1' = (-b_1 y_1 - k_1 (x_1 - L_1) + k_2 (x_2 - x_1 - L_2))/m_1$

$x_2' = y_2$

$y_2' = (-b_2 y_2 - k_2 (x_2 - x_1 - L_2))/m_2$

现在这些方程已经可以被我们用 Python 实现。

以下代码定义了方程组的“右侧”（也称为向量场）。我选择将定义向量场的函数放在它自己的模块中（即在它自己的文件中），但这不是必需的。请注意，函数的参数被配置为与函数 odeint 一起使用：时间 t 是第二个参数。

在 [1]
def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled spring-mass system.

    Arguments:
        w :  vector of the state variables:
                  w = [x1,y1,x2,y2]
        t :  time
        p :  vector of the parameters:
                  p = [m1,m2,k1,k2,L1,L2,b1,b2]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, L1, L2, b1, b2 = p

    # Create f = (x1',y1',x2',y2'):
    f = [y1,
         (-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
         y2,
         (-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
    return f

接下来，这里有一个脚本使用 odeint 来求解给定参数值、初始条件和时间间隔的方程。该脚本将点写入文件 'two_springs.dat'。

在 [2]
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint

# 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.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
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(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

with open('two_springs.dat', 'w') as f:
    # Print & save the solution.
    for t1, w1 in zip(t, wsol):
        print >> f, t1, w1[0], w1[1], w1[2], w1[3]

以下脚本使用 Matplotlib 来绘制由上述脚本生成的解。

在 [3]
# Plot the solution that was generated

from numpy import loadtxt
from pylab import figure, plot, xlabel, grid, hold, legend, title, savefig
from matplotlib.font_manager import FontProperties

t, x1, xy, x2, y2 = loadtxt('two_springs.dat', unpack=True)

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

xlabel('t')
grid(True)
hold(True)
lw = 1

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

legend((r'$x_1$', r'$x_2$'), prop=FontProperties(size=16))
title('Mass Displacements for the\nCoupled Spring-Mass System')
savefig('two_springs.png', dpi=100)

章节作者：WarrenWeckesser，Warren Weckesser

附件