我是ouint函数，但我不确定函数是否允许这样做 . 我已经从Scipy Cookbook修改了这个例子，以便耦合的两个体系可能会发生碰撞 .
图片链接
图碰撞
想象一下
具体而言， spring 常数和质量已经变得更弱更轻 . 您还会注意到函数内部使用的'if'语句试图限制'm1'的行程，而不是'x2'的初始输入 . 对不起，我没有足够的声望点来发布图表，但你会清楚地看到两个质量可以在解决方案的某些部分占用相同的空间 .
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
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]
if y1 > x2:
y1 == x2
else:
y1 == y1
return f
# Parameter values
# Masses:
m1 = 0.5
m2 = 1.5
# Spring constants
k1 = 0.1
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 = 4.25
y2 = 0.0
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 5.0
numpoints = 2500
# 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)
plt.plot(t, wsol[:, 0], 'b', label='theta(t)')
plt.plot(t, wsol[:, 2], 'g', label='omega(t)')
plt.show()
我不确定我的变量是否设置不正确，或者odeint函数是否不能接受我试图添加的这些类型的“限制” . 理想情况下，如果可以设置“上边界”，如果质量的位移超过该值，则可以将其设置为等于最大排序限制 .
提前感谢您的帮助 .