如何在Rust的ode_solvers中提取积分全程状态用于绘图？

1. 获取积分过程的中间值与最终值 #

Dop853步进器完成integrate()后，可通过两个核心方法提取数据：

• x_out()：返回所有记录的时间点切片（类型为&[Time]）
• y_out()：返回对应每个时间点的状态向量切片（类型为&[State]）

修改你的main函数匹配块，添加数据提取逻辑：

match res {
    Ok(stats) => {
        println!("stats: {:?}", stats);
        
        // 提取时间序列与状态数据
        let times = stepper.x_out();
        let states = stepper.y_out();
        
        // 示例：打印初始与最终状态
        println!("初始状态: {:?}", states.first().unwrap());
        println!("最终状态: {:?}", states.last().unwrap());
        
        // 提取位置分量（x/y/z对应状态前三个元素）
        let positions: Vec<(f64, f64, f64)> = states.iter()
            .map(|s| (s[0], s[1], s[2]))
            .collect();
    },
    Err(_) => println!("An error occurred."),
}

2. 绘制轨迹示例（搭配plotters库） #

要生成轨迹图，可使用plotters绘图库。首先在Cargo.toml添加依赖：

[dependencies]
ode_solvers = "0.16.0"
plotters = "0.3.5"

完整可运行代码如下：

use ode_solvers::{self, dop853::*};
use plotters::prelude::*;

type State = ode_solvers::Vector6<f64>;
type Time = f64;

struct KeplerOrbit {
    mu: f64,
}

impl ode_solvers::System<Time, State> for KeplerOrbit {
    fn system(&self, _t: Time, y: &State, dy: &mut State) {
        let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt();

        dy[0] = y[3];
        dy[1] = y[4];
        dy[2] = y[5];
        dy[3] = -self.mu * y[0] / r.powi(3);
        dy[4] = -self.mu * y[1] / r.powi(3);
        dy[5] = -self.mu * y[2] / r.powi(3);
    }
}

fn main() {
    let system = KeplerOrbit { mu: 398600.435436 };

    let a: f64 = 20000.0;
    let period = 2.0 * std::f64::consts::PI * (a.powi(3) / system.mu).sqrt();
    let y0 = State::new(
        -5007.248417988539,
        -1444.918140151374,
        3628.534606178356,
        0.717716656891,
        -10.224093784269,
        0.748229399696,
    );

    // 第四个参数为记录步长，控制中间点的密度，值越小保存的点越多
    let mut stepper = Dop853::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10);
    let res = stepper.integrate();

    match res {
        Ok(stats) => {
            println!("stats: {:?}", stats);

            let states = stepper.y_out();
            // 提取X-Y平面的位置数据用于绘图
            let xy_positions: Vec<(f64, f64)> = states.iter().map(|s| (s[0], s[1])).collect();

            // 创建绘图后端并生成图片
            let root = BitMapBackend::new("kepler_orbit.png", (800, 600)).into_drawing_area();
            root.fill(&WHITE).unwrap();
            let mut chart = ChartBuilder::on(&root)
                .caption("Kepler Orbit (X-Y Plane)", ("sans-serif", 20).into_font())
                .x_label_area_size(40)
                .y_label_area_size(40)
                .build_cartesian_2d(-25000.0..25000.0, -25000.0..25000.0)
                .unwrap();

            chart.configure_mesh().draw().unwrap();

            // 绘制轨迹线
            chart
                .draw_series(LineSeries::new(xy_positions, &BLUE))
                .unwrap()
                .label("Orbit Trajectory")
                .legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &BLUE));

            chart.configure_series_labels().draw().unwrap();

            root.present().unwrap();
            println!("轨迹图已保存为 kepler_orbit.png");
        }
        Err(_) => println!("An error occurred."),
    }
}

关键注意事项 #

• Dop853::new的第四个参数是记录步长，决定每隔多少时间保存一次状态数据，若需要更精细的轨迹，可减小该值。
• y_out()中的每个State是6维向量：前三个元素为位置(x/y/z)，后三个为速度(vx/vy/vz)，可根据需求提取对应分量。
• 运行代码后会在当前目录生成kepler_orbit.png，展示X-Y平面的开普勒轨道轨迹。

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