如何从Runge Kutta 4模拟的抛射点数据中提取最高点曲线？

下面提供两种实现思路，一种是通过脚本预处理数据提取最高点，另一种是直接用gnuplot完成数据处理与绘图。

一、预处理数据：用脚本提取所有轨迹的最高点 #

如果你的轨迹数据是离散采样的，每个轨迹的最高点可近似为该轨迹中y值最大的点（采样密度足够时误差可忽略）。根据数据存储方式选择对应脚本：

单轨迹文件（每个发射角度对应单独文件） #

假设轨迹文件命名为traj_theta_0.dat、traj_theta_1.dat等，用Python脚本遍历提取最高点：

import os
import numpy as np

output_file = "max_points.dat"
with open(output_file, 'w') as f_out:
    for filename in os.listdir('./trajectories'):
        if filename.startswith('traj_theta_') and filename.endswith('.dat'):
            data = np.loadtxt(os.path.join('./trajectories', filename))
            x, y = data[:, 0], data[:, 1]
            max_idx = np.argmax(y)
            f_out.write(f"{x[max_idx]} {y[max_idx]}\n")

之后用gnuplot绘制：

plot "max_points.dat" with linespoints lw 3 lc rgb "black" title "最高点曲线", \
     "traj_*.dat" with lines lc rgb "light-gray" title "各轨迹"

合并轨迹文件（所有角度数据在同一文件，格式为x y theta） #

先按角度分组，每组提取y最大值对应的点：

import numpy as np

data = np.loadtxt('all_trajectories.dat')
thetas = np.unique(data[:, 2])
with open('max_points.dat', 'w') as f_out:
    for theta in thetas:
        traj = data[data[:,2] == theta]
        max_idx = np.argmax(traj[:,1])
        f_out.write(f"{traj[max_idx,0]} {traj[max_idx,1]}\n")

二、直接用gnuplot处理（无需额外脚本） #

gnuplot内置的stats命令和循环功能可直接从原始数据中提取最高点：

单轨迹文件场景 #

# 初始化存储数组（根据轨迹数量调整数组大小）
array max_x[100]
array max_y[100]
count = 0

# 遍历所有轨迹文件
do for [file in system("ls traj_theta_*.dat")] {
    stats file using 1:2 nooutput
    count = count + 1
    max_x[count] = STATS_pos_max_y_x  # y最大值对应的x
    max_y[count] = STATS_max_y        # y最大值
}

# 绘制所有轨迹与最高点曲线
plot for [file in system("ls traj_theta_*.dat")] file with lines lc rgb "light-gray" notitle, \
     '-' with lines lw 3 lc rgb "black" title "最高点曲线"
do for [i=1:count] {
    print max_x[i], max_y[i]
}
e

合并轨迹文件场景（数据按角度连续分组） #

# 获取所有唯一发射角度
stats "all_trajectories.dat" using 3 nooutput
array thetas[STATS_blocks]
do for [i=1:STATS_blocks] {
    thetas[i] = STATS_block_value[i]
}

# 收集每个角度的最高点
array max_x[STATS_blocks]
array max_y[STATS_blocks]
do for [i=1:STATS_blocks] {
    stats "all_trajectories.dat" using 1:2 every ::i::i nooutput
    max_x[i] = STATS_pos_max_y_x
    max_y[i] = STATS_max_y
}

# 绘图
plot "all_trajectories.dat" using 1:2 with lines lc rgb "light-gray" notitle, \
     '-' with lines lw 3 lc rgb "black" title "最高点曲线"
do for [i=1:STATS_blocks] {
    print max_x[i], max_y[i]
}
e

注意事项 #

• 若采样稀疏，直接取y最大值会有误差，建议在RK4计算时，当y方向速度接近0时增加采样密度，或用三次样条插值拟合轨迹后再找精确最高点。
• 用gnuplot处理时，确保数据格式规范，比如单轨迹文件命名有规律，合并文件中同一角度的轨迹数据连续。

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