如何串联3个简谐振子？基于Python简谐振子代码实现级联

修改逻辑 #

原代码实现的是单个无驱动自由简谐振子的数值求解，要实现3个相同振子串联，核心调整如下：

• 单个振子包含「位移、速度」2个状态变量，3个串联振子共需维护6个状态变量，顺序为[x0, v0, x1, v1, x2, v2]，分别对应第1/2/3个振子的位移和速度
• 第一个振子保持原自由振荡的动力学规则，初始条件沿用原设置：初始位移0，初始速度2；其余两个振子初始状态为静止在平衡位置（位移0、速度0）
• 串联振子为耦合动力学系统：后一个振子的受力由相邻前一个振子的位移差决定，所有振子保持相同的omega_sq参数
• 一次性求解整个耦合系统的状态演化，最终同时绘制3个振子的位移随时间变化曲线方便对比

修改后完整代码 #

import numpy as np 
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

# 时间轴参数和原代码保持一致
t = np.linspace(0, 15, 1000)
omega_sq = 1
# 初始状态：[第1个振子位移, 第1个振子速度, 第2个振子位移, 第2个振子速度, 第3个振子位移, 第3个振子速度]
init_state = [0, 2, 0, 0, 0, 0]

def series_oscillator(t, state):
    x0, v0, x1, v1, x2, v2 = state
    # 位移的导数是速度，速度的导数是加速度（由弹簧弹力提供）
    dx0 = v0
    dv0 = -omega_sq * x0  # 第一个振子一端固定，回复力仅和自身位移相关
    dx1 = v1
    dv1 = -omega_sq * (x1 - x0)  # 第二个振子受力由前后位移差决定
    dx2 = v2
    dv2 = -omega_sq * (x2 - x1)  # 第三个振子受力由前后位移差决定
    return [dx0, dv0, dx1, dv1, dx2, dv2]

# 求解耦合系统
res = solve_ivp(series_oscillator, t_span=[t.min(), t.max()], y0=init_state, t_eval=t)

# 绘制三个振子的位移曲线
plt.plot(t, res.y[0], label='Oscillator 1')
plt.plot(t, res.y[2], label='Oscillator 2')
plt.plot(t, res.y[4], label='Oscillator 3')
plt.ylabel("Position")
plt.xlabel("Time")
plt.title('3 Series Coupled Simple Harmonic Oscillators', fontsize=16)
plt.legend()
plt.show()

补充说明 #

如果你需要的是信号级联模式（即前一个振子的输出作为独立外驱动信号输入给后级振子，而非物理弹簧刚性耦合），只需要把后级振子的动力学方程改为受迫振动形式，分三次依次求解即可：先算出第一个振子的位移序列，通过插值把该序列作为驱动项输入第二个振子求解，再用第二个振子的输出驱动第三个振子。

内容的提问来源于stack exchange，提问作者j123