用Python实现弹簧振子运动仿真四阶龙格-库塔法实战指南在工程仿真和物理系统建模中弹簧振子是最基础却最具代表性的动力学模型之一。从机械减震器到分子动力学理解如何准确模拟弹簧振子的运动规律是每个工程师和物理学家的必修课。本文将带你用Python实现一个完整的弹簧振子仿真系统重点介绍四阶龙格-库塔法这一经典数值解法在实际物理问题中的应用。1. 弹簧振子系统的数学建模1.1 物理系统分析考虑一个典型的弹簧-质量-阻尼系统它由以下要素组成质量为m的物体弹性系数为k的弹簧阻尼系数为c的阻尼器当系统偏离平衡位置时会受到三种力的作用弹簧的恢复力$F_{spring} -kx$阻尼器的阻力$F_{damping} -cv$可能的外部驱动力$F_{ext}(t)$根据牛顿第二定律系统的运动方程可以表示为$$ m\frac{d^2x}{dt^2} c\frac{dx}{dt} kx F_{ext}(t) $$1.2 微分方程的标准化处理为了便于数值求解我们需要将这个二阶常微分方程转化为一阶方程组。设$$ \begin{cases} y_1 x \ y_2 \frac{dx}{dt} \end{cases} $$则原方程可以转化为$$ \begin{cases} \frac{dy_1}{dt} y_2 \ \frac{dy_2}{dt} \frac{1}{m}(F_{ext}(t) - cy_2 - ky_1) \end{cases} $$这种形式非常适合用龙格-库塔法进行数值求解。2. 四阶龙格-库塔法原理详解2.1 方法的基本思想四阶龙格-库塔法(RK4)是一种单步显式数值积分方法通过计算四个不同位置的斜率估计值然后进行加权平均来获得更高精度的解。对于一般的一阶微分方程$$ \frac{dy}{dt} f(t, y) $$RK4的迭代公式为$$ y_{n1} y_n \frac{h}{6}(k_1 2k_2 2k_3 k_4) $$其中$k_1 f(t_n, y_n)$$k_2 f(t_n h/2, y_n hk_1/2)$$k_3 f(t_n h/2, y_n hk_2/2)$$k_4 f(t_n h, y_n hk_3)$2.2 应用于弹簧振子系统对于我们的弹簧振子系统需要同时处理两个一阶方程。定义$$ \mathbf{y} \begin{bmatrix} y_1 \ y_2 \end{bmatrix}, \quad \mathbf{f}(t, \mathbf{y}) \begin{bmatrix} y_2 \ \frac{1}{m}(F_{ext}(t) - cy_2 - ky_1) \end{bmatrix} $$则RK4的向量形式为$$ \mathbf{y}_{n1} \mathbf{y}_n \frac{h}{6}(\mathbf{k}_1 2\mathbf{k}_2 2\mathbf{k}_3 \mathbf{k}_4) $$其中各$\mathbf{k}$的计算方式与标量情况类似。3. Python实现完整代码3.1 基础实现import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class SpringMassDamper: def __init__(self, m, c, k): self.m m # 质量(kg) self.c c # 阻尼系数(N·s/m) self.k k # 弹簧系数(N/m) def equation(self, t, y, F_ext0): 定义系统微分方程 y1, y2 y dy1dt y2 dy2dt (F_ext - self.c*y2 - self.k*y1) / self.m return np.array([dy1dt, dy2dt]) def rk4_step(self, t, y, h, F_ext0): 单步RK4求解 k1 self.equation(t, y, F_ext) k2 self.equation(t h/2, y h/2*k1, F_ext) k3 self.equation(t h/2, y h/2*k2, F_ext) k4 self.equation(t h, y h*k3, F_ext) y_new y h/6 * (k1 2*k2 2*k3 k4) return y_new def solve(self, t_span, y0, h, F_ext_funcNone): 求解整个时间区间的运动 t_start, t_end t_span num_steps int((t_end - t_start) / h) t_values np.linspace(t_start, t_end, num_steps1) y_values np.zeros((num_steps1, 2)) y_values[0] y0 for i in range(num_steps): t t_values[i] F_ext F_ext_func(t) if F_ext_func else 0 y_values[i1] self.rk4_step(t, y_values[i], h, F_ext) return t_values, y_values3.2 可视化与结果分析def plot_results(t, y, title): 绘制位移和速度随时间变化的曲线 plt.figure(figsize(12, 6)) plt.subplot(2, 1, 1) plt.plot(t, y[:, 0], b-, labelDisplacement (m)) plt.ylabel(Displacement (m)) plt.title(title) plt.grid(True) plt.legend() plt.subplot(2, 1, 2) plt.plot(t, y[:, 1], r-, labelVelocity (m/s)) plt.xlabel(Time (s)) plt.ylabel(Velocity (m/s)) plt.grid(True) plt.legend() plt.tight_layout() plt.show() def animate_motion(t, y): 创建动画展示质量块的运动 fig, ax plt.subplots(figsize(10, 6)) ax.set_xlim(-2, 2) ax.set_ylim(-0.5, 0.5) ax.grid(True) spring, ax.plot([], [], b-, lw2) mass plt.Rectangle((0, -0.2), 0.4, 0.4, fcr) ax.add_patch(mass) time_text ax.text(0.02, 0.95, , transformax.transAxes) def init(): spring.set_data([], []) mass.set_xy((0, -0.2)) time_text.set_text() return spring, mass, time_text def update(frame): x y[frame, 0] # 绘制弹簧(简化表示) spring_x np.linspace(-1.5, x-0.2, 20) spring_y 0.1 * np.sin(10*(spring_x1.5)/(x1.3)) spring.set_data(spring_x, spring_y) mass.set_xy((x, -0.2)) time_text.set_text(fTime {t[frame]:.2f}s, Displacement {x:.3f}m) return spring, mass, time_text ani FuncAnimation(fig, update, frameslen(t), init_funcinit, blitTrue, interval50) plt.close() return ani4. 实际应用案例4.1 自由振动分析考虑一个无阻尼自由振动系统质量m 1 kg弹簧系数k 10 N/m初始位移x0 0.5 m初始速度v0 0 m/s# 创建系统实例 system SpringMassDamper(m1, c0, k10) # 初始条件 y0 np.array([0.5, 0]) # [初始位移, 初始速度] # 求解 t, y system.solve(t_span(0, 10), y0y0, h0.01) # 绘制结果 plot_results(t, y, Free Vibration (No Damping)) # 生成动画 ani animate_motion(t, y) from IPython.display import HTML HTML(ani.to_jshtml())4.2 阻尼振动分析增加阻尼系数c 0.5 N·s/m# 创建阻尼系统实例 damped_system SpringMassDamper(m1, c0.5, k10) # 求解 t_damped, y_damped damped_system.solve(t_span(0, 10), y0y0, h0.01) # 比较自由振动和阻尼振动 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelNo damping) plt.plot(t_damped, y_damped[:, 0], r-, labelWith damping (c0.5)) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Free and Damped Vibration) plt.legend() plt.grid(True) plt.show()4.3 受迫振动分析考虑一个周期性外力作用下的系统def external_force(t): # 幅值为2N频率为1Hz的正弦外力 return 2 * np.sin(2*np.pi*1*t) # 创建系统实例(带阻尼) forced_system SpringMassDamper(m1, c0.2, k10) # 求解 t_forced, y_forced forced_system.solve( t_span(0, 10), y0y0, h0.01, F_ext_funcexternal_force ) # 绘制结果 plot_results(t_forced, y_forced, Forced Vibration with Damping) # 绘制外力与位移响应 plt.figure(figsize(10, 5)) plt.plot(t_forced, external_force(t_forced), g-, labelExternal Force (N)) plt.plot(t_forced, y_forced[:, 0], b-, labelDisplacement (m)) plt.xlabel(Time (s)) plt.ylabel(Amplitude) plt.title(Forced Vibration: Force Input vs System Response) plt.legend() plt.grid(True) plt.show()5. 方法验证与误差分析5.1 解析解对比对于无阻尼自由振动系统我们可以得到解析解$$ x(t) A\cos(\omega_n t) B\sin(\omega_n t) $$其中$\omega_n \sqrt{k/m}$是固有频率。# 计算解析解 omega_n np.sqrt(10) # sqrt(k/m) A y0[0] # 初始位移 B y0[1]/omega_n # 初始速度/omega_n analytical_solution A * np.cos(omega_n * t) B * np.sin(omega_n * t) # 比较数值解和解析解 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelNumerical (RK4)) plt.plot(t, analytical_solution, r--, labelAnalytical) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Numerical and Analytical Solutions) plt.legend() plt.grid(True) plt.show() # 计算误差 error np.abs(y[:, 0] - analytical_solution) plt.figure(figsize(10, 5)) plt.plot(t, error, m-) plt.xlabel(Time (s)) plt.ylabel(Absolute Error (m)) plt.title(Numerical Error Over Time) plt.grid(True) plt.show()5.2 步长对精度的影响# 测试不同步长下的误差 step_sizes [0.1, 0.05, 0.01, 0.005] max_errors [] for h in step_sizes: t_h, y_h system.solve(t_span(0, 10), y0y0, hh) analytical_h A * np.cos(omega_n * t_h) B * np.sin(omega_n * t_h) max_error np.max(np.abs(y_h[:, 0] - analytical_h)) max_errors.append(max_error) # 绘制误差随步长变化 plt.figure(figsize(10, 5)) plt.loglog(step_sizes, max_errors, bo-) plt.xlabel(Step Size (s)) plt.ylabel(Maximum Absolute Error (m)) plt.title(Error vs Step Size for RK4 Method) plt.grid(True, whichboth) plt.show()6. 高级应用与扩展6.1 参数敏感性分析研究不同参数对系统响应的影响# 测试不同阻尼系数的影响 damping_coeffs [0, 0.2, 0.5, 1, 2] plt.figure(figsize(10, 6)) for c in damping_coeffs: test_system SpringMassDamper(m1, cc, k10) t_test, y_test test_system.solve(t_span(0, 10), y0y0, h0.01) plt.plot(t_test, y_test[:, 0], labelfc{c} N·s/m) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(System Response with Different Damping Coefficients) plt.legend() plt.grid(True) plt.show()6.2 非线性弹簧系统考虑弹簧力与位移不成正比的非线性情况class NonlinearSpringSystem(SpringMassDamper): def equation(self, t, y, F_ext0): y1, y2 y dy1dt y2 # 非线性弹簧力F_spring -k*y1 - α*y1^3 alpha 2 # 非线性系数 dy2dt (F_ext - self.c*y2 - self.k*y1 - alpha*y1**3) / self.m return np.array([dy1dt, dy2dt]) # 创建非线性系统实例 nonlinear_system NonlinearSpringSystem(m1, c0.1, k10) # 求解 t_nonlinear, y_nonlinear nonlinear_system.solve(t_span(0, 10), y0y0, h0.01) # 比较线性和非线性响应 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelLinear system) plt.plot(t_nonlinear, y_nonlinear[:, 0], g-, labelNonlinear system) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Linear and Nonlinear System Responses) plt.legend() plt.grid(True) plt.show()6.3 多自由度系统扩展将方法扩展到耦合的弹簧-质量系统class CoupledSpringMassSystem: def __init__(self, m1, m2, k1, k2, c1, c2): self.m1 m1 self.m2 m2 self.k1 k1 self.k2 k2 self.c1 c1 self.c2 c2 def equation(self, t, y): x1, v1, x2, v2 y # 两个质量块的加速度 a1 (-self.k1*x1 - self.c1*v1 self.k2*(x2-x1) self.c2*(v2-v1)) / self.m1 a2 (-self.k2*(x2-x1) - self.c2*(v2-v1)) / self.m2 return np.array([v1, a1, v2, a2]) def rk4_step(self, t, y, h): k1 self.equation(t, y) k2 self.equation(t h/2, y h/2*k1) k3 self.equation(t h/2, y h/2*k2) k4 self.equation(t h, y h*k3) y_new y h/6 * (k1 2*k2 2*k3 k4) return y_new def solve(self, t_span, y0, h): t_start, t_end t_span num_steps int((t_end - t_start) / h) t_values np.linspace(t_start, t_end, num_steps1) y_values np.zeros((num_steps1, 4)) y_values[0] y0 for i in range(num_steps): y_values[i1] self.rk4_step(t_values[i], y_values[i], h) return t_values, y_values # 创建耦合系统实例 coupled_system CoupledSpringMassSystem( m11, m21.5, k110, k28, c10.1, c20.05 ) # 初始条件 [x1, v1, x2, v2] y0_coupled np.array([0.5, 0, -0.3, 0]) # 求解 t_coupled, y_coupled coupled_system.solve(t_span(0, 15), y0y0_coupled, h0.01) # 绘制结果 plt.figure(figsize(12, 6)) plt.plot(t_coupled, y_coupled[:, 0], b-, labelMass 1 displacement) plt.plot(t_coupled, y_coupled[:, 2], r-, labelMass 2 displacement) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Coupled Spring-Mass System Response) plt.legend() plt.grid(True) plt.show()
用Python搞定物理仿真：四阶龙格-库塔法求解弹簧振子运动方程（附完整代码）
用Python实现弹簧振子运动仿真四阶龙格-库塔法实战指南在工程仿真和物理系统建模中弹簧振子是最基础却最具代表性的动力学模型之一。从机械减震器到分子动力学理解如何准确模拟弹簧振子的运动规律是每个工程师和物理学家的必修课。本文将带你用Python实现一个完整的弹簧振子仿真系统重点介绍四阶龙格-库塔法这一经典数值解法在实际物理问题中的应用。1. 弹簧振子系统的数学建模1.1 物理系统分析考虑一个典型的弹簧-质量-阻尼系统它由以下要素组成质量为m的物体弹性系数为k的弹簧阻尼系数为c的阻尼器当系统偏离平衡位置时会受到三种力的作用弹簧的恢复力$F_{spring} -kx$阻尼器的阻力$F_{damping} -cv$可能的外部驱动力$F_{ext}(t)$根据牛顿第二定律系统的运动方程可以表示为$$ m\frac{d^2x}{dt^2} c\frac{dx}{dt} kx F_{ext}(t) $$1.2 微分方程的标准化处理为了便于数值求解我们需要将这个二阶常微分方程转化为一阶方程组。设$$ \begin{cases} y_1 x \ y_2 \frac{dx}{dt} \end{cases} $$则原方程可以转化为$$ \begin{cases} \frac{dy_1}{dt} y_2 \ \frac{dy_2}{dt} \frac{1}{m}(F_{ext}(t) - cy_2 - ky_1) \end{cases} $$这种形式非常适合用龙格-库塔法进行数值求解。2. 四阶龙格-库塔法原理详解2.1 方法的基本思想四阶龙格-库塔法(RK4)是一种单步显式数值积分方法通过计算四个不同位置的斜率估计值然后进行加权平均来获得更高精度的解。对于一般的一阶微分方程$$ \frac{dy}{dt} f(t, y) $$RK4的迭代公式为$$ y_{n1} y_n \frac{h}{6}(k_1 2k_2 2k_3 k_4) $$其中$k_1 f(t_n, y_n)$$k_2 f(t_n h/2, y_n hk_1/2)$$k_3 f(t_n h/2, y_n hk_2/2)$$k_4 f(t_n h, y_n hk_3)$2.2 应用于弹簧振子系统对于我们的弹簧振子系统需要同时处理两个一阶方程。定义$$ \mathbf{y} \begin{bmatrix} y_1 \ y_2 \end{bmatrix}, \quad \mathbf{f}(t, \mathbf{y}) \begin{bmatrix} y_2 \ \frac{1}{m}(F_{ext}(t) - cy_2 - ky_1) \end{bmatrix} $$则RK4的向量形式为$$ \mathbf{y}_{n1} \mathbf{y}_n \frac{h}{6}(\mathbf{k}_1 2\mathbf{k}_2 2\mathbf{k}_3 \mathbf{k}_4) $$其中各$\mathbf{k}$的计算方式与标量情况类似。3. Python实现完整代码3.1 基础实现import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class SpringMassDamper: def __init__(self, m, c, k): self.m m # 质量(kg) self.c c # 阻尼系数(N·s/m) self.k k # 弹簧系数(N/m) def equation(self, t, y, F_ext0): 定义系统微分方程 y1, y2 y dy1dt y2 dy2dt (F_ext - self.c*y2 - self.k*y1) / self.m return np.array([dy1dt, dy2dt]) def rk4_step(self, t, y, h, F_ext0): 单步RK4求解 k1 self.equation(t, y, F_ext) k2 self.equation(t h/2, y h/2*k1, F_ext) k3 self.equation(t h/2, y h/2*k2, F_ext) k4 self.equation(t h, y h*k3, F_ext) y_new y h/6 * (k1 2*k2 2*k3 k4) return y_new def solve(self, t_span, y0, h, F_ext_funcNone): 求解整个时间区间的运动 t_start, t_end t_span num_steps int((t_end - t_start) / h) t_values np.linspace(t_start, t_end, num_steps1) y_values np.zeros((num_steps1, 2)) y_values[0] y0 for i in range(num_steps): t t_values[i] F_ext F_ext_func(t) if F_ext_func else 0 y_values[i1] self.rk4_step(t, y_values[i], h, F_ext) return t_values, y_values3.2 可视化与结果分析def plot_results(t, y, title): 绘制位移和速度随时间变化的曲线 plt.figure(figsize(12, 6)) plt.subplot(2, 1, 1) plt.plot(t, y[:, 0], b-, labelDisplacement (m)) plt.ylabel(Displacement (m)) plt.title(title) plt.grid(True) plt.legend() plt.subplot(2, 1, 2) plt.plot(t, y[:, 1], r-, labelVelocity (m/s)) plt.xlabel(Time (s)) plt.ylabel(Velocity (m/s)) plt.grid(True) plt.legend() plt.tight_layout() plt.show() def animate_motion(t, y): 创建动画展示质量块的运动 fig, ax plt.subplots(figsize(10, 6)) ax.set_xlim(-2, 2) ax.set_ylim(-0.5, 0.5) ax.grid(True) spring, ax.plot([], [], b-, lw2) mass plt.Rectangle((0, -0.2), 0.4, 0.4, fcr) ax.add_patch(mass) time_text ax.text(0.02, 0.95, , transformax.transAxes) def init(): spring.set_data([], []) mass.set_xy((0, -0.2)) time_text.set_text() return spring, mass, time_text def update(frame): x y[frame, 0] # 绘制弹簧(简化表示) spring_x np.linspace(-1.5, x-0.2, 20) spring_y 0.1 * np.sin(10*(spring_x1.5)/(x1.3)) spring.set_data(spring_x, spring_y) mass.set_xy((x, -0.2)) time_text.set_text(fTime {t[frame]:.2f}s, Displacement {x:.3f}m) return spring, mass, time_text ani FuncAnimation(fig, update, frameslen(t), init_funcinit, blitTrue, interval50) plt.close() return ani4. 实际应用案例4.1 自由振动分析考虑一个无阻尼自由振动系统质量m 1 kg弹簧系数k 10 N/m初始位移x0 0.5 m初始速度v0 0 m/s# 创建系统实例 system SpringMassDamper(m1, c0, k10) # 初始条件 y0 np.array([0.5, 0]) # [初始位移, 初始速度] # 求解 t, y system.solve(t_span(0, 10), y0y0, h0.01) # 绘制结果 plot_results(t, y, Free Vibration (No Damping)) # 生成动画 ani animate_motion(t, y) from IPython.display import HTML HTML(ani.to_jshtml())4.2 阻尼振动分析增加阻尼系数c 0.5 N·s/m# 创建阻尼系统实例 damped_system SpringMassDamper(m1, c0.5, k10) # 求解 t_damped, y_damped damped_system.solve(t_span(0, 10), y0y0, h0.01) # 比较自由振动和阻尼振动 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelNo damping) plt.plot(t_damped, y_damped[:, 0], r-, labelWith damping (c0.5)) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Free and Damped Vibration) plt.legend() plt.grid(True) plt.show()4.3 受迫振动分析考虑一个周期性外力作用下的系统def external_force(t): # 幅值为2N频率为1Hz的正弦外力 return 2 * np.sin(2*np.pi*1*t) # 创建系统实例(带阻尼) forced_system SpringMassDamper(m1, c0.2, k10) # 求解 t_forced, y_forced forced_system.solve( t_span(0, 10), y0y0, h0.01, F_ext_funcexternal_force ) # 绘制结果 plot_results(t_forced, y_forced, Forced Vibration with Damping) # 绘制外力与位移响应 plt.figure(figsize(10, 5)) plt.plot(t_forced, external_force(t_forced), g-, labelExternal Force (N)) plt.plot(t_forced, y_forced[:, 0], b-, labelDisplacement (m)) plt.xlabel(Time (s)) plt.ylabel(Amplitude) plt.title(Forced Vibration: Force Input vs System Response) plt.legend() plt.grid(True) plt.show()5. 方法验证与误差分析5.1 解析解对比对于无阻尼自由振动系统我们可以得到解析解$$ x(t) A\cos(\omega_n t) B\sin(\omega_n t) $$其中$\omega_n \sqrt{k/m}$是固有频率。# 计算解析解 omega_n np.sqrt(10) # sqrt(k/m) A y0[0] # 初始位移 B y0[1]/omega_n # 初始速度/omega_n analytical_solution A * np.cos(omega_n * t) B * np.sin(omega_n * t) # 比较数值解和解析解 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelNumerical (RK4)) plt.plot(t, analytical_solution, r--, labelAnalytical) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Numerical and Analytical Solutions) plt.legend() plt.grid(True) plt.show() # 计算误差 error np.abs(y[:, 0] - analytical_solution) plt.figure(figsize(10, 5)) plt.plot(t, error, m-) plt.xlabel(Time (s)) plt.ylabel(Absolute Error (m)) plt.title(Numerical Error Over Time) plt.grid(True) plt.show()5.2 步长对精度的影响# 测试不同步长下的误差 step_sizes [0.1, 0.05, 0.01, 0.005] max_errors [] for h in step_sizes: t_h, y_h system.solve(t_span(0, 10), y0y0, hh) analytical_h A * np.cos(omega_n * t_h) B * np.sin(omega_n * t_h) max_error np.max(np.abs(y_h[:, 0] - analytical_h)) max_errors.append(max_error) # 绘制误差随步长变化 plt.figure(figsize(10, 5)) plt.loglog(step_sizes, max_errors, bo-) plt.xlabel(Step Size (s)) plt.ylabel(Maximum Absolute Error (m)) plt.title(Error vs Step Size for RK4 Method) plt.grid(True, whichboth) plt.show()6. 高级应用与扩展6.1 参数敏感性分析研究不同参数对系统响应的影响# 测试不同阻尼系数的影响 damping_coeffs [0, 0.2, 0.5, 1, 2] plt.figure(figsize(10, 6)) for c in damping_coeffs: test_system SpringMassDamper(m1, cc, k10) t_test, y_test test_system.solve(t_span(0, 10), y0y0, h0.01) plt.plot(t_test, y_test[:, 0], labelfc{c} N·s/m) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(System Response with Different Damping Coefficients) plt.legend() plt.grid(True) plt.show()6.2 非线性弹簧系统考虑弹簧力与位移不成正比的非线性情况class NonlinearSpringSystem(SpringMassDamper): def equation(self, t, y, F_ext0): y1, y2 y dy1dt y2 # 非线性弹簧力F_spring -k*y1 - α*y1^3 alpha 2 # 非线性系数 dy2dt (F_ext - self.c*y2 - self.k*y1 - alpha*y1**3) / self.m return np.array([dy1dt, dy2dt]) # 创建非线性系统实例 nonlinear_system NonlinearSpringSystem(m1, c0.1, k10) # 求解 t_nonlinear, y_nonlinear nonlinear_system.solve(t_span(0, 10), y0y0, h0.01) # 比较线性和非线性响应 plt.figure(figsize(10, 5)) plt.plot(t, y[:, 0], b-, labelLinear system) plt.plot(t_nonlinear, y_nonlinear[:, 0], g-, labelNonlinear system) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Comparison of Linear and Nonlinear System Responses) plt.legend() plt.grid(True) plt.show()6.3 多自由度系统扩展将方法扩展到耦合的弹簧-质量系统class CoupledSpringMassSystem: def __init__(self, m1, m2, k1, k2, c1, c2): self.m1 m1 self.m2 m2 self.k1 k1 self.k2 k2 self.c1 c1 self.c2 c2 def equation(self, t, y): x1, v1, x2, v2 y # 两个质量块的加速度 a1 (-self.k1*x1 - self.c1*v1 self.k2*(x2-x1) self.c2*(v2-v1)) / self.m1 a2 (-self.k2*(x2-x1) - self.c2*(v2-v1)) / self.m2 return np.array([v1, a1, v2, a2]) def rk4_step(self, t, y, h): k1 self.equation(t, y) k2 self.equation(t h/2, y h/2*k1) k3 self.equation(t h/2, y h/2*k2) k4 self.equation(t h, y h*k3) y_new y h/6 * (k1 2*k2 2*k3 k4) return y_new def solve(self, t_span, y0, h): t_start, t_end t_span num_steps int((t_end - t_start) / h) t_values np.linspace(t_start, t_end, num_steps1) y_values np.zeros((num_steps1, 4)) y_values[0] y0 for i in range(num_steps): y_values[i1] self.rk4_step(t_values[i], y_values[i], h) return t_values, y_values # 创建耦合系统实例 coupled_system CoupledSpringMassSystem( m11, m21.5, k110, k28, c10.1, c20.05 ) # 初始条件 [x1, v1, x2, v2] y0_coupled np.array([0.5, 0, -0.3, 0]) # 求解 t_coupled, y_coupled coupled_system.solve(t_span(0, 15), y0y0_coupled, h0.01) # 绘制结果 plt.figure(figsize(12, 6)) plt.plot(t_coupled, y_coupled[:, 0], b-, labelMass 1 displacement) plt.plot(t_coupled, y_coupled[:, 2], r-, labelMass 2 displacement) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(Coupled Spring-Mass System Response) plt.legend() plt.grid(True) plt.show()
