阻尼振动的可视化 包括源码和推导

flyfish

牛顿第二定律（加速度定律）
胡克定律（Hooke‘s Law）

阻尼振动是指在振动系统中，由于阻力或能量损耗导致振动幅度随时间减小的现象。

左边为无阻尼，右边为有阻尼，有阻尼的摆动会逐渐停止。

对于简单的阻尼振动系统，其运动方程为：
$$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+kx=0$$
在摆动系统中，上述方程的对应形式为：
$$\frac{d\omega }{dt}=-\frac{b}{l}\omega -\frac{g}{l}\sin (\theta )$$

1. 简单摆动方程（无阻尼） ：$$\frac{d\theta }{dt}=\omega$$

$$\frac{d\omega }{dt}=-\frac{g}{l}\sin (\theta )$$

其中：

$g$（重力加速度）对应代码中的 GRAVITY = 9.8。

$l$（摆长）对应代码中的 LENGTH = 1.0。

2. 阻尼摆动方程 ：$$\frac{d\theta }{dt}=\omega$$

$$\frac{d\omega }{dt}=-\frac{b}{l}\omega -\frac{g}{l}\sin (\theta )$$

其中：

$g$（重力加速度）对应代码中的 GRAVITY = 9.8。
$l$（摆长）对应代码中的 LENGTH = 1.0。
$b$（阻尼系数）对应代码中的 DAMPING_COEFFICIENT = 0.2。

阻尼振动系统的推导

1. 牛顿第二定律

对于一个质量 $m$ 的物体，其受力情况可以用牛顿第二定律来描述：
$F=m\frac{d^{2}x}{d{t}^2}$
其中 $x(t)$ 是物体的位置函数， $\frac{d^{2}x}{d{t}^2}$ 是加速度。

2. 力的组成

在一个简单的阻尼振动系统中，力由以下几部分组成：

弹力 ：遵循胡克定律，弹力与位移成正比，并且方向与位移相反。
$F_{\text{spring}}=-kx$
其中 $k$ 是弹簧常数， $x$ 是位移。

阻尼力 ：阻尼力与速度成正比，并且方向与速度相反。
$F_{\text{damping}}=-b\frac{dx}{dt}$
其中 $b$ 是阻尼系数， $\frac{dx}{dt}$ 是速度。

重力和其他外力 ：对于水平振动的简单系统，通常可以忽略重力和其他外力。如果考虑垂直振动，重力可以被包含在弹力的静态部分中。

3. 总受力

总受力为上述各力的叠加：
$F=F_{\text{spring}}+F_{\text{damping}}=-kx-b\frac{dx}{dt}$

4. 代入牛顿第二定律

根据牛顿第二定律：
$m\frac{d^{2}x}{d{t}^2}=-kx-b\frac{dx}{dt}$

5. 标准形式

将上述方程整理为标准形式：
$m\frac{d^{2}x}{d{t}^2}+b\frac{dx}{dt}+kx=0$这个方程描述了一个质量 $m$ 的物体在弹簧常数 $k$ 和阻尼系数 $b$ 下的阻尼振动。可以将其简化为无量纲形式，定义以下参数：
固有角频率（不含阻尼）： ${\omega }_0=\sqrt{\frac{k}{m}}$

阻尼比： $\gamma =\frac{b}{2m}$
于是方程可以写成：
$\frac{d^{2}x}{d{t}^2}+2\gamma \frac{dx}{dt}+{\omega }_0^{2}x=0$

运动方程

对于阻尼振动系统，方程描述的是阻尼摆动：
$\frac{d^{2}x}{d{t}^2}+2\gamma \frac{dx}{dt}+{\omega }_0^{2}x=0$

动画代码中的公式

之前在动画代码中使用的阻尼摆动方程是：
$\frac{d\theta }{dt}=\omega$
$\frac{d\omega }{dt}=-b\omega -\frac{g}{l}\sin (\theta )$这里，$\theta$ 是摆角，$\omega$ 是角速度：
$b$ 是阻尼系数，对应代码中的 DAMPING_COEFFICIENT。

$\frac{g}{l}$ 是摆动系统中的恢复力系数，对应代码中的 GRAVITY / LENGTH。

可视化源码实现

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import odeint
from math import sin, cos

# Constants
GRAVITY = 9.8
LENGTH = 1.0
DAMPING_COEFFICIENT = 0.2

def pendulum_equations_no_damping(state, t, length):
    """Equations for a simple pendulum without damping."""
    theta, omega = state
    dtheta_dt = omega
    domega_dt = -GRAVITY / length * sin(theta)
    return dtheta_dt, domega_dt

def pendulum_equations_damping(state, t, length, damping):
    """Equations for a simple pendulum with damping."""
    theta, omega = state
    dtheta_dt = omega
    domega_dt = -damping * omega - GRAVITY / length * sin(theta)
    return dtheta_dt, domega_dt

def solve_pendulum(equations, initial_state, t, args):
    """Solves the pendulum ODE."""
    return odeint(equations, initial_state, t, args=args)

def calculate_trajectory(track, length):
    """Calculates the x and y coordinates of the pendulum bob."""
    x_data = length * np.sin(track[:, 0])
    y_data = -length * np.cos(track[:, 0])
    return x_data, y_data

def initialize_animation():
    """Initializes the animation plot."""
    for ax in axes:
        ax.set_xlim(-1.5 * LENGTH, 1.5 * LENGTH)
        ax.set_ylim(-1.5 * LENGTH, 1.5 * LENGTH)
        ax.grid()
    time_text_no_damping.set_text('')
    time_text_damping.set_text('')
    return lines + time_texts

def update_animation(frame):
    """Updates the animation for each frame."""
    lines[0].set_data([0, x_data_no_damping[frame]], [0, y_data_no_damping[frame]])
    lines[1].set_data([0, x_data_damping[frame]], [0, y_data_damping[frame]])
    time_text_no_damping.set_text(f'time = {frame * 0.1:.1f}s')
    time_text_damping.set_text(f'time = {frame * 0.1:.1f}s')
    return lines + time_texts

# Time array
t = np.arange(0, 20, 0.1)

# Initial state (theta, omega)
initial_state = (1.0, 0)

# Solve ODE for pendulum with and without damping
track_no_damping = solve_pendulum(pendulum_equations_no_damping, initial_state, t, args=(LENGTH,))
track_damping = solve_pendulum(pendulum_equations_damping, initial_state, t, args=(LENGTH, DAMPING_COEFFICIENT))

# Calculate trajectory
x_data_no_damping, y_data_no_damping = calculate_trajectory(track_no_damping, LENGTH)
x_data_damping, y_data_damping = calculate_trajectory(track_damping, LENGTH)

# Create plot
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
lines = [axes[0].plot([], [], 'o-', lw=2, color='blue')[0], 
         axes[1].plot([], [], 'o-', lw=2, color='green')[0]]
time_template = 'time = %.1fs'
time_text_no_damping = axes[0].text(0.05, 0.9, '', transform=axes[0].transAxes)
time_text_damping = axes[1].text(0.05, 0.9, '', transform=axes[1].transAxes)
time_texts = [time_text_no_damping, time_text_damping]

# Add titles
axes[0].set_title('Pendulum without Damping')
axes[1].set_title('Pendulum with Damping')

# Create animation
ani = animation.FuncAnimation(
    fig, update_animation, frames=len(x_data_no_damping), init_func=initialize_animation, interval=50
)

# Save animation
ani.save('double_pendulum.gif', writer='imagemagick', fps=100)
plt.show()