Python，调用scipy.integrate.solve_ivp有条件的二度spring-质量系统

Question

我正在尝试使用 scipy.integrate.

中的 solve_ivp 函数来解决像这样的 2 自由度 spring 质量系统

图像的link如下。这是我的第一个post。所以堆栈溢出不允许我 post 图片。

这里，当质量 m1 在负区域超过 d 时， K1 变成两倍，即当

-infinity < x1 <= -1, k1 = 14
-1 < x1 <= infinity, k1 = 7

这是我为 m1

编写的没有第二个 spring 的代码

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def F(t, y):
##########################################################################################################################################################################################
    m1 = 3
    m2 = 5

    k1 = 7
    k2 = 9

    c1 = 1
    c2 = 2

    f1 = 40*np.cos(3*t)
    f2 = 0
##########################################################################################################################################################################################

    M = np.array([
        [m1, 0],
        [0, m2]
    ])

    C = np.array([
        [c1 + c2, -c2],
        [-c2, c2]
    ])

    K = np.array([
        [k1 + k2, -k2],
        [-k2, k2]
    ])

    A = np.vstack([
        np.hstack([
            np.zeros((2, 2)), 
            np.eye(2, 2)
        ]), 
        np.hstack([
            -np.matmul(np.linalg.inv(M), K), 
            -np.matmul(np.linalg.inv(M), C)
        ])
    ])

    B = np.vstack([
        np.zeros((2, 2)), 
        np.linalg.inv(M)
    ])

    F = np.array([
        [f1],
        [f2]
    ])

    yvec = np.array([[y[i] for i in range(4)]]).T

    ydot = np.matmul(A, yvec) + np.matmul(B, F)

    return ydot

####################################################################################################################################################################################
start_time = 0
end_time = 60
delta_t = 0.01

initial_position_m_1 = 2
initial_velocity_m_1 = 1

initial_position_m_2 = 5
initial_velocity_m_2 = 3
####################################################################################################################################################################################

TE = np.arange(start_time, end_time, delta_t)

time_interval = np.array([start_time, end_time])
initial_conditions = np.array([initial_position_m_1, initial_position_m_2, initial_velocity_m_1, initial_velocity_m_2])

sol = solve_ivp(F, time_interval, initial_conditions, t_eval = TE, vectorized=True, method = 'RK45')


T = sol.t
Y = sol.y
dis_m_1 = sol.y[0]
dis_m_2 = sol.y[1]
vel_m_1 = sol.y[2]
vel_m_2 = sol.y[3]

plt.plot(T, dis_m_1)
plt.plot(T, np.full((len(dis_m_1)), -1, dtype='int32'))
plt.title("Displacement of m1 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, dis_m_2)
plt.title("Displacement of m2 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, vel_m_1)
plt.title("Velocity of m1 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, vel_m_2)
plt.title("Velocity of m2 vs time")
plt.xlabel("Time")
plt.show()

这是m1没有第二个spring的位移图。

link of the image

只要质量 m1 低于水平线，k1 就应该等于 14。如何在以下代码的 return 语句中为 m1

应用 if else 条件 2 springs

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def F(t, y):
##########################################################################################################################################################################################
    m1 = 3
    m2 = 5

    k1 = 7
    k2 = 9

    c1 = 1
    c2 = 2

    f1 = 40*np.cos(3*t)
    f2 = 0
##########################################################################################################################################################################################

    M = np.array([
        [m1, 0],
        [0, m2]
    ])

    C = np.array([
        [c1 + c2, -c2],
        [-c2, c2]
    ])

    K = np.array([
        [k1 + k2, -k2],
        [-k2, k2]
    ])

    KK = np.array([
        [2*k1 + k2, -k2],
        [-k2, k2]
    ])

    A = np.vstack([
        np.hstack([
            np.zeros((2, 2)), 
            np.eye(2, 2)
        ]), 
        np.hstack([
            -np.matmul(np.linalg.inv(M), K), 
            -np.matmul(np.linalg.inv(M), C)
        ])
    ])

    AA = np.vstack([
        np.hstack([
            np.zeros((2, 2)), 
            np.eye(2, 2)
        ]), 
        np.hstack([
            -np.matmul(np.linalg.inv(M), KK), 
            -np.matmul(np.linalg.inv(M), C)
        ])
    ])


    B = np.vstack([
        np.zeros((2, 2)), 
        np.linalg.inv(M)
    ])

    F = np.array([
        [f1],
        [f2]
    ])

    yvec = np.array([[y[i] for i in range(4)]]).T

    ydot1 = np.matmul(A, yvec) + np.matmul(B, F)    # k1 = 7 
    ydot2 = np.matmul(AA, yvec) + np.matmul(B, F)   # k1 = 14

    return ???????

####################################################################################################################################################################################
start_time = 0
end_time = 60
delta_t = 0.01

initial_position_m_1 = 2
initial_velocity_m_1 = 1

initial_position_m_2 = 5
initial_velocity_m_2 = 3
####################################################################################################################################################################################

TE = np.arange(start_time, end_time, delta_t)

time_interval = np.array([start_time, end_time])
initial_conditions = np.array([initial_position_m_1, initial_position_m_2, initial_velocity_m_1, initial_velocity_m_2])

sol = solve_ivp(F, time_interval, initial_conditions, t_eval = TE, vectorized=True, method = 'RK45')


T = sol.t
Y = sol.y
dis_m_1 = sol.y[0]
dis_m_2 = sol.y[1]
vel_m_1 = sol.y[2]
vel_m_2 = sol.y[3]

plt.plot(T, dis_m_1)
plt.plot(T, np.full((len(dis_m_1)), -1, dtype='int32'))
plt.title("Displacement of m1 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, dis_m_2)
plt.title("Displacement of m2 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, vel_m_1)
plt.title("Velocity of m1 vs time")
plt.xlabel("Time")
plt.show()

plt.plot(T, vel_m_2)
plt.title("Velocity of m2 vs time")
plt.xlabel("Time")
plt.show()

我觉得应该是这样的。

if displacement_of_m1 <= -1:
    return ydot1
else:
    return ydot2

有人可以帮忙吗？谢谢。

Answer 1

应该是

if y[0] < -1:
    return ydot2
return ydot1

此处 y[0]、y[1] ... 最多 y[n] 存在，它可能会根据 vectorized=True 发生变化。有时它会是 [x1, x2, x3, ..., v1, v2, v3...] 或 [x1, v1, x2, v2, x3, v3, ...]

Python，调用scipy.integrate.solve_ivp有条件的二度spring-质量系统

Python, calling scipy.integrate.solve_ivp with conditions for a second degree spring-mass system

python

scipy

ode

solver

differential-equations