diff --git a/tools/NumericalMethods.py b/tools/NumericalMethods.py new file mode 100644 index 0000000..3528611 --- /dev/null +++ b/tools/NumericalMethods.py @@ -0,0 +1,331 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation + +# 设置中文字体 +plt.rcParams['font.sans-serif'] = ['SimHei'] +plt.rcParams['axes.unicode_minus'] = False + +# 简谐振子系统参数 +k = 1.0 # 劲度系数 +m = 1.0 # 质量 +omega = np.sqrt(k/m) # 角频率 + +def potential_energy(x): + """势能: (1/2)*k*x^2""" + return 0.5 * k * x**2 + +def kinetic_energy(v): + """动能: (1/2)*m*v^2""" + return 0.5 * m * v**2 + +def total_energy(x, v): + """总能量""" + return potential_energy(x) + kinetic_energy(v) + +def acceleration(x): + """加速度: a = -(k/m)*x""" + return -omega**2 * x + +# ============ 四种数值方法 ============ + +def euler_explicit(x0, v0, dt, n_steps): + """欧拉显式 (前向欧拉)""" + x = np.zeros(n_steps) + v = np.zeros(n_steps) + t = np.zeros(n_steps) + x[0] = x0 + v[0] = v0 + t[0] = 0 + for i in range(n_steps - 1): + x[i+1] = x[i] + v[i] * dt + v[i+1] = v[i] + acceleration(x[i]) * dt + t[i+1] = t[i] + dt + return t, x, v + + +def euler_implicit(x0, v0, dt, n_steps): + """欧拉隐式 (后向欧拉)""" + x = np.zeros(n_steps) + v = np.zeros(n_steps) + t = np.zeros(n_steps) + x[0] = x0 + v[0] = v0 + t[0] = 0 + for i in range(n_steps - 1): + # 解析解：对简谐振子联立方程可直接得出 + x[i+1] = (x[i] + v[i] * dt) / (1 + omega**2 * dt**2) + v[i+1] = v[i] + acceleration(x[i+1]) * dt + t[i+1] = t[i] + dt + return t, x, v + + +def midpoint_method(x0, v0, dt, n_steps): + """中点差分法 (二阶 Runge-Kutta / 改进欧拉)""" + x = np.zeros(n_steps) + v = np.zeros(n_steps) + t = np.zeros(n_steps) + x[0] = x0 + v[0] = v0 + t[0] = 0 + for i in range(n_steps - 1): + x_mid = x[i] + v[i] * dt / 2 + v_mid = v[i] + acceleration(x[i]) * dt / 2 + x[i+1] = x[i] + v_mid * dt + v[i+1] = v[i] + acceleration(x_mid) * dt + t[i+1] = t[i] + dt + return t, x, v + + +def leapfrog_method(x0, v0, dt, n_steps): + """蛙跳法 (Leapfrog / Störmer-Verlet) + + 核心思想：位置和速度交错半步更新 + v_{n+1/2} = v_{n-1/2} + a(x_n) * dt <- 速度在半整数时刻 + x_{n+1} = x_n + v_{n+1/2} * dt <- 位置在整数时刻 + + 启动：用欧拉半步得到 v_{1/2}，之后进入交错循环。 + 输出时将速度还原到整数时刻：v_n = (v_{n-1/2} + v_{n+1/2}) / 2 + """ + x = np.zeros(n_steps) + v = np.zeros(n_steps) # 整数时刻速度（仅用于输出/能量计算） + t = np.zeros(n_steps) + + x[0] = x0 + v[0] = v0 + t[0] = 0 + + # 用欧拉半步初始化 v_{1/2} + v_half = v0 + acceleration(x0) * dt / 2 + + for i in range(n_steps - 1): + # 更新位置（整数时刻） + x[i+1] = x[i] + v_half * dt + t[i+1] = t[i] + dt + + # 更新半步速度（用新位置的加速度） + v_half_new = v_half + acceleration(x[i+1]) * dt + + # 整数时刻速度（用于能量计算和输出）：取两个半步的平均 + v[i+1] = (v_half + v_half_new) / 2 + + v_half = v_half_new + + return t, x, v + + +def exact_solution(x0, v0, t): + """精确解""" + A = np.sqrt(x0**2 + (v0 / omega)**2) + phi = np.arctan2(-v0 / omega, x0) + x_exact = A * np.cos(omega * t + phi) + v_exact = -A * omega * np.sin(omega * t + phi) + return x_exact, v_exact + + +# ============ 参数设置 ============ +dt = 0.05 # 时间步长 +T_total = 10.0 # 总时间 +n_steps = int(T_total / dt) + 1 + +# 三个初始条件 (x0, v0) +initial_conditions = [ + (1.0, 0.0), # 从最大位移释放 + (0.0, 1.0), # 从平衡位置给初速度 + (0.8, 0.6) # 混合初始条件 +] + +# 四种方法：(名称, 函数, 线型颜色, alpha) +methods = [ + ('欧拉显式', euler_explicit, 'r-', 0.8), + ('欧拉隐式', euler_implicit, 'b-', 0.8), + ('中点差分法', midpoint_method, 'g-', 0.8), + ('蛙跳法', leapfrog_method, 'm-', 0.8), +] + +colors_ic = ['#e63946', '#457b9d', '#2d6a4f'] # 三种初始条件的颜色 +labels_ic = ['从最大位移释放', '平衡位置初速度', '混合条件'] + +# 四种方法对应的绘图颜色（用于位移/速度时序图） +method_colors = ['r', 'b', 'g', 'm'] + +print(f"时间步长 dt = {dt}s, 步数 = {n_steps}") + +# ============ 计算结果 ============ +results = {} +for method_name, method_func, color_style, alpha in methods: + results[method_name] = {} + for idx, (x0, v0) in enumerate(initial_conditions): + t_arr, x_arr, v_arr = method_func(x0, v0, dt, n_steps) + E = total_energy(x_arr, v_arr) + results[method_name][idx] = { + 't': t_arr, 'x': x_arr, 'v': v_arr, 'E': E, + 'E_initial': E[0], 'E_final': E[-1], + 'E_drift': (E[-1] - E[0]) / E[0] * 100 + } + +# 精确解（高密度时间轴） +exact_t = np.linspace(0, T_total, 1000) + +# ============ 图形布局 ============ +fig = plt.figure(figsize=(18, 12)) +fig.suptitle('简谐振子数值方法四方比较', fontsize=16, fontweight='bold') + +gs = fig.add_gridspec(3, 4, hspace=0.40, wspace=0.30) +ax_phase = fig.add_subplot(gs[0, :]) # 第1行全宽：相空间 +ax_energy = fig.add_subplot(gs[1, :]) # 第2行全宽：能量演化 +ax_traj = fig.add_subplot(gs[2, 0]) # 第3行：位置-时间 +ax_velocity = fig.add_subplot(gs[2, 1]) # 第3行：速度-时间 +ax_err_bar = fig.add_subplot(gs[2, 2]) # 第3行：能量误差柱状图 +ax_animation= fig.add_subplot(gs[2, 3]) # 第3行：相空间动画 + +# ============ 1. 相空间轨迹 ============ +ax_phase.set_title('相空间轨迹比较 (x vs v)', fontsize=13) +ax_phase.set_xlabel('位置 x', fontsize=11) +ax_phase.set_ylabel('速度 v', fontsize=11) +ax_phase.grid(True, alpha=0.3) +ax_phase.axhline(y=0, color='k', linestyle='-', alpha=0.3) +ax_phase.axvline(x=0, color='k', linestyle='-', alpha=0.3) + +line_styles_ic = ['-', '--', ':'] +method_line_colors = ['r', 'b', 'g', 'm'] + +for midx, (x0, v0) in enumerate(initial_conditions): + x_exact, v_exact = exact_solution(x0, v0, exact_t) + ax_phase.plot(x_exact, v_exact, 'k-', alpha=0.25, linewidth=1, + label='精确解' if midx == 0 else "") + for mi, (method_name, _, _, _) in enumerate(methods): + res = results[method_name][midx] + ax_phase.plot(res['x'], res['v'], + linestyle=['-', '--', ':', '-.'][midx], + color=method_line_colors[mi], alpha=0.75, linewidth=1.4, + label=method_name if midx == 0 else "") + +handles, labels_leg = ax_phase.get_legend_handles_labels() +by_label = dict(zip(labels_leg, handles)) +ax_phase.legend(by_label.values(), by_label.keys(), + loc='upper right', fontsize=8, ncol=2) + +# ============ 2. 能量演化（归一化） ============ +ax_energy.set_title('归一化能量随时间演化 E(t)/E(0)', fontsize=13) +ax_energy.set_xlabel('时间 t (s)', fontsize=11) +ax_energy.set_ylabel('E / E0', fontsize=11) +ax_energy.grid(True, alpha=0.3) + +energy_linestyles = ['-', '--', ':', '-.'] +for mi, (method_name, _, _, _) in enumerate(methods): + res = results[method_name][0] # 只用第一个初始条件展示能量 + E_norm = res['E'] / res['E_initial'] + ax_energy.plot(res['t'], E_norm, + linestyle=energy_linestyles[mi], + color=method_line_colors[mi], + linewidth=1.8, alpha=0.9, label=method_name) + +ax_energy.axhline(y=1.0, color='k', linestyle='-', alpha=0.5, linewidth=0.8) +ax_energy.set_ylim(0.85, 1.70) +ax_energy.legend(loc='upper left', fontsize=9) + +# ============ 3. 位置-时间 (第一初始条件) ============ +idx_show = 0 +ax_traj.set_title(f'位置-时间 ({labels_ic[idx_show]})', fontsize=11) +ax_traj.set_xlabel('时间 t (s)', fontsize=10) +ax_traj.set_ylabel('位置 x', fontsize=10) +ax_traj.grid(True, alpha=0.3) + +x_ex, v_ex = exact_solution(*initial_conditions[idx_show], exact_t) +ax_traj.plot(exact_t, x_ex, 'k-', alpha=0.5, linewidth=1.2, label='精确解') +for mi, (method_name, _, _, _) in enumerate(methods): + res = results[method_name][idx_show] + ax_traj.plot(res['t'], res['x'], + color=method_line_colors[mi], linewidth=1.4, + alpha=0.85, label=method_name) +ax_traj.legend(loc='upper right', fontsize=7) + +# ============ 4. 速度-时间 (第一初始条件) ============ +ax_velocity.set_title(f'速度-时间 ({labels_ic[idx_show]})', fontsize=11) +ax_velocity.set_xlabel('时间 t (s)', fontsize=10) +ax_velocity.set_ylabel('速度 v', fontsize=10) +ax_velocity.grid(True, alpha=0.3) + +ax_velocity.plot(exact_t, v_ex, 'k-', alpha=0.5, linewidth=1.2, label='精确解') +for mi, (method_name, _, _, _) in enumerate(methods): + res = results[method_name][idx_show] + ax_velocity.plot(res['t'], res['v'], + color=method_line_colors[mi], linewidth=1.4, + alpha=0.85, label=method_name) +ax_velocity.legend(loc='upper right', fontsize=7) + +# ============ 5. 能量误差柱状图 ============ +ax_err_bar.set_title('能量漂移百分比 |dE/E0|*100%', fontsize=11) +ax_err_bar.set_ylabel('漂移 (%)', fontsize=10) +ax_err_bar.grid(True, alpha=0.3, axis='y') + +method_names_short = ['显式', '隐式', '中点', '蛙跳'] +bar_colors = method_line_colors +drifts = [abs(results[mn][0]['E_drift']) for mn, *_ in methods] +bars = ax_err_bar.bar(method_names_short, drifts, color=bar_colors, alpha=0.8, width=0.5) +for bar, val in zip(bars, drifts): + ax_err_bar.text(bar.get_x() + bar.get_width() / 2, + bar.get_height() + 0.3, + f'{val:.2f}%', ha='center', va='bottom', fontsize=9) +ax_err_bar.set_ylim(0, max(drifts) * 1.2 + 2) + +# ============ 6. 动画：四种方法相空间对比 (同一初始条件) ============ +ax_animation.set_xlim(-1.5, 1.5) +ax_animation.set_ylim(-1.5, 1.5) +ax_animation.set_title('四种方法同步释放 (相空间)', fontsize=11) +ax_animation.set_xlabel('位置 x', fontsize=10) +ax_animation.set_ylabel('速度 v', fontsize=10) +ax_animation.grid(True, alpha=0.3) +ax_animation.axhline(y=0, color='k', linestyle='-', alpha=0.3) +ax_animation.axvline(x=0, color='k', linestyle='-', alpha=0.3) + +anim_points = [] +anim_trails = [] +anim_idx = 0 # 选第一组初始条件，比较四种方法的同步演化 +for mi, (method_name, _, _, _) in enumerate(methods): + pt, = ax_animation.plot([], [], 'o', color=method_line_colors[mi], + markersize=7, label=method_name) + tr, = ax_animation.plot([], [], '-', color=method_line_colors[mi], alpha=0.35, linewidth=1) + anim_points.append(pt) + anim_trails.append(tr) + +x_anim_exact, v_anim_exact = exact_solution(*initial_conditions[anim_idx], exact_t) +ax_animation.plot(x_anim_exact, v_anim_exact, 'k-', alpha=0.2, linewidth=0.8, label='精确解') + +ax_animation.legend(loc='upper right', fontsize=8) + +def animate(frame): + for mi, (method_name, _, _, _) in enumerate(methods): + res = results[method_name][anim_idx] + if frame < len(res['x']): + anim_points[mi].set_data([res['x'][frame]], [res['v'][frame]]) + start = max(0, frame - 50) + anim_trails[mi].set_data(res['x'][start:frame+1], res['v'][start:frame+1]) + return anim_points + anim_trails + +anim = FuncAnimation(fig, animate, frames=n_steps, interval=20, blit=True) + +# ============ 打印能量误差统计 ============ +print("\n" + "=" * 65) +print("能量守恒误差分析 (相对漂移百分比，初始条件：从最大位移释放)") +print("=" * 65) +for method_name, _, _, _ in methods: + res0 = results[method_name][0] + print(f" {method_name:10s}：{res0['E_drift']:+.4f}%") + +print("\n" + "=" * 65) +print("所有初始条件下的能量漂移") +print("=" * 65) +for method_name, _, _, _ in methods: + print(f"\n{method_name}:") + for idx, (x0, v0) in enumerate(initial_conditions): + res = results[method_name][idx] + print(f" {labels_ic[idx]}: {res['E_drift']:+.4f}%") + +# 显示图形 +plt.tight_layout() +plt.show() + +# 可选：保存动画 +# anim.save('harmonic_oscillator_4methods.gif', writer='pillow', fps=30)