我正在用 Rust 实现 Runge-Kutta-Fehlberg 5(4) 方法，所以我必须跟踪我试图建模的系统的当前状态，根据当前状态执行几项计算，并根据这些计算的结果更新状态。
在此过程中，我遇到了所有权问题，因为nalgebra我使用的类型没有实现copy。编译器建议我使用clone数据。我最终.clone()在一个循环中写了四条语句，我一方面想知道这会对性能造成什么影响，另一方面想知道是否有其他可行的替代方案。
您可以在下面找到相关功能。
pub fn rkf54(
initial_time: f64,
initial_state: State<f64>,
initial_step_size: f64,
control: bool,
tolerance: f64,
final_time: f64,
f: &dyn Fn(f64, &State<f64>) -> State<f64>)
-> Vec<State<f64>> {
let mut output: Vec<State<f64>> = vec![];
let mut time = initial_time;
let mut state = initial_state;
let mut step_size = initial_step_size;
while time < final_time {
if control {
// CLONE 1
step_size = rkf_control(time, state.clone(), initial_step_size, tolerance, f);
}
// CLONE 2
let stages = rkf_stages(time, &state.clone(), step_size, A, C, f);
// CLONE 3
state = rkf_step(&state.clone(), step_size, &stages, B_PRIMARY);
time += step_size;
// CLONE 4
output.push(state.clone())
}
let stages = rkf_stages(final_time, &state.clone(), step_size, A, C, f);
state = rkf_step(&state.clone(), step_size, &stages, B_PRIMARY);
output.push(state.clone());
return output;
}
重复克隆会明显影响性能吗？如果我理解正确的话，每次迭代都会在内存的不同区域复制相同的数据4次。
如果是这样，我该如何避免这种克隆/分配？
这是一个完整的工作示例以及类型和常量定义。
#![allow(dead_code)]
extern crate nalgebra; // 0.33.2
use nalgebra::DVector; // 0.33.2
use std::f64::INFINITY;
type Stages<T, const D: usize> = [DVector<T>; D];
pub type State<T> = DVector<T>;
pub type StageCoefficients<T, const D: usize> = [&'static [T]; D];
pub type Weights<T, const D: usize> = [T; D];
pub const A: StageCoefficients<f64, 6> = [
&[],
&[1./4.],
&[3./32., 9./32.],
&[1932./2197., -7200./2197., 7296./2197.],
&[439./216., 8., 3680./513., -845./4104.],
&[-8./27., 2., -3544./2565., 1859./4104., -11./40.]
];
pub const C: Weights<f64, 6> = [0.0, 1./4., 3./8., 12./13., 1., 1./2. ];
pub const B_PRIMARY: Weights<f64, 6> = [16./135., 0., 6656./12825., 28561./56430., -9./50., 2./55. ];
pub const B_EMBEDDED: Weights<f64, 6> = [25./216., 0., 1408./2565., 2197./4104., -1./5., 0. ];
const K: f64 = 0.02;
const B0: f64 = 3.5;
const B1: f64 = 16.72;
pub fn f(t:f64, state: &State<f64>) -> State<f64> {
let x = state[0];
let y = state [1];
let x_dot = y;
let y_dot = (-K*y) - (x.powf(3.)) + (B0) + (B1 * t.cos());
return State::from(vec![x_dot, y_dot]);
}
pub fn rkf54(
initial_time: f64,
initial_state: State<f64>,
initial_step_size: f64,
control: bool,
tolerance: f64,
final_time: f64,
f: &dyn Fn(f64, &State<f64>) -> State<f64>)
-> Vec<State<f64>> {
let mut output: Vec<State<f64>> = vec![];
let mut time = initial_time;
let mut state = initial_state;
let mut step_size = initial_step_size;
while time < final_time {
if control {
// CLONE 1
step_size = rkf_control(time, state.clone(), initial_step_size, tolerance, f);
}
// CLONE 2
let stages = rkf_stages(time, &state.clone(), step_size, A, C, f);
// CLONE 3
state = rkf_step(&state.clone(), step_size, &stages, B_PRIMARY);
time += step_size;
// CLONE 4
output.push(state.clone())
}
let stages = rkf_stages(final_time, &state.clone(), step_size, A, C, f);
state = rkf_step(&state.clone(), step_size, &stages, B_PRIMARY);
output.push(state.clone());
return output;
}
fn rkf_stages<const D: usize>(
time: f64,
state: &State<f64>,
step_size: f64,
stage_coefficients: StageCoefficients<f64, D>,
time_fractions: Weights<f64, D>,
f: &dyn Fn(f64, &State<f64>) -> State<f64>
) -> Stages<f64, D> {
let mut k: Stages<f64, D> = std::array::from_fn(|_| DVector::zeros(state.nrows()));
for i in 0..stage_coefficients.len() {
let mut x: State<f64> = DVector::zeros(state.nrows());
for j in 0..stage_coefficients[i].len() {
x += &k[j] * stage_coefficients[i][j];
}
x = state + step_size*(x);
let t = time + (step_size*time_fractions[i]);
k[i] = f(t,&x);
}
return k;
}
fn rkf_step<const D: usize>(state: &State<f64>, step_size: f64, stages: &Stages<f64, D>, weights: Weights<f64, D>) -> State<f64> {
let mut x_out = 0.0 * DVector::zeros(state.nrows());
for i in 0..A.len() {
x_out += &stages[i] * weights[i];
}
x_out *= step_size;
x_out += state;
return x_out;
}
fn rkf_control(time: f64, state: State<f64>, initial_step_size: f64, tolerance: f64,
f: &dyn Fn(f64, &State<f64>) -> State<f64>) -> f64 {
let mut error_estimate = INFINITY;
let mut step_size = initial_step_size;
while error_estimate > tolerance {
let k = rkf_stages(time, &state, step_size, A, C, f);
let x_primary = rkf_step(&state, step_size, &k, B_PRIMARY);
let x_embedded = rkf_step(&state, step_size, &k, B_EMBEDDED);
error_estimate = (x_primary - x_embedded).norm();
step_size = 0.9 * step_size * (tolerance / error_estimate).powf(1. / 5.);
}
return step_size;
}
fn main() {
let initial_time = 0.;
let initial_state = State::from(vec![2., 2.]);
let initial_step_size = 1.;
let control = true;
let tolerance = 10e-3;
let final_time = 1.;
let output = rkf54(initial_time, initial_state, initial_step_size, control, tolerance, final_time, &f);
for state in output {
println!("{state:1.3e}");
}
}
