pub fn rk4_step(y: &[f64], t: f64, dt: f64, f: &dyn Fn(f64, &[f64]) -> Vec<f64>) -> Vec<f64> {
let n = y.len();
let k1 = f(t, y);
let mut y2 = vec![0.0; n];
for i in 0..n {
y2[i] = y[i] + 0.5 * dt * k1[i];
}
let k2 = f(t + 0.5 * dt, &y2);
let mut y3 = vec![0.0; n];
for i in 0..n {
y3[i] = y[i] + 0.5 * dt * k2[i];
}
let k3 = f(t + 0.5 * dt, &y3);
let mut y4 = vec![0.0; n];
for i in 0..n {
y4[i] = y[i] + dt * k3[i];
}
let k4 = f(t + dt, &y4);
let mut result = vec![0.0; n];
for i in 0..n {
result[i] = y[i] + dt / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}
result
}
pub fn rk4_integrate(
y0: &[f64],
t_start: f64,
t_end: f64,
n_steps: usize,
f: &dyn Fn(f64, &[f64]) -> Vec<f64>,
) -> Vec<Vec<f64>> {
let dt = (t_end - t_start) / n_steps as f64;
let mut trajectory = Vec::with_capacity(n_steps + 1);
let mut y = y0.to_vec();
trajectory.push(y.clone());
let mut t = t_start;
for _ in 0..n_steps {
y = rk4_step(&y, t, dt, f);
t += dt;
trajectory.push(y.clone());
}
trajectory
}
pub fn leapfrog_step(
positions: &mut [f64],
velocities: &mut [f64],
dt: f64,
acceleration: &dyn Fn(&[f64]) -> Vec<f64>,
) {
let n = positions.len();
let a_initial = acceleration(positions);
for i in 0..n {
velocities[i] += 0.5 * dt * a_initial[i];
}
for i in 0..n {
positions[i] += dt * velocities[i];
}
let a_final = acceleration(positions);
for i in 0..n {
velocities[i] += 0.5 * dt * a_final[i];
}
}
pub fn leapfrog_integrate(
positions: &mut [f64],
velocities: &mut [f64],
dt: f64,
n_steps: usize,
acceleration: &dyn Fn(&[f64]) -> Vec<f64>,
) -> Vec<(Vec<f64>, Vec<f64>)> {
let mut trajectory = Vec::with_capacity(n_steps + 1);
trajectory.push((positions.to_vec(), velocities.to_vec()));
for _ in 0..n_steps {
leapfrog_step(positions, velocities, dt, acceleration);
trajectory.push((positions.to_vec(), velocities.to_vec()));
}
trajectory
}
pub fn root_bisection(
f: &dyn Fn(f64) -> f64,
mut a: f64,
mut b: f64,
tol: f64,
max_iter: usize,
) -> f64 {
for _ in 0..max_iter {
let mid = 0.5 * (a + b);
if (b - a).abs() < tol {
return mid;
}
if f(a) * f(mid) <= 0.0 {
b = mid;
} else {
a = mid;
}
}
0.5 * (a + b)
}
pub fn root_brent(
f: &dyn Fn(f64) -> f64,
mut a: f64,
mut b: f64,
tol: f64,
max_iter: usize,
) -> f64 {
let mut fa = f(a);
let mut fb = f(b);
if fa * fb > 0.0 {
return 0.5 * (a + b);
}
let mut c = a;
let mut fc = fa;
let mut d = b - a;
let mut e = d;
for _ in 0..max_iter {
if fb * fc > 0.0 {
c = a;
fc = fa;
d = b - a;
e = d;
}
if fc.abs() < fb.abs() {
a = b;
b = c;
c = a;
fa = fb;
fb = fc;
fc = fa;
}
let tol1 = 2.0 * f64::EPSILON * b.abs() + 0.5 * tol;
let m = 0.5 * (c - b);
if m.abs() <= tol1 || fb.abs() < 1e-50 {
return b;
}
if e.abs() >= tol1 && fa.abs() > fb.abs() {
let s_val = fb / fa;
let (p, q_val) = if (a - c).abs() < 1e-50 {
(2.0 * m * s_val, 1.0 - s_val)
} else {
let q_inner = fa / fc;
let r = fb / fc;
(
s_val * (2.0 * m * q_inner * (q_inner - r) - (b - a) * (r - 1.0)),
(q_inner - 1.0) * (r - 1.0) * (s_val - 1.0),
)
};
let q_val = if p > 0.0 { -q_val } else { q_val };
let p = p.abs();
if 2.0 * p < (3.0 * m * q_val - (tol1 * q_val).abs()).min(e * q_val) {
e = d;
d = p / q_val;
} else {
d = m;
e = m;
}
} else {
d = m;
e = m;
}
a = b;
fa = fb;
if d.abs() > tol1 {
b += d;
} else {
b += if m > 0.0 { tol1 } else { -tol1 };
}
fb = f(b);
}
b
}
pub fn root_newton(
f: &dyn Fn(f64) -> f64,
df: &dyn Fn(f64) -> f64,
mut x: f64,
tol: f64,
max_iter: usize,
) -> f64 {
for _ in 0..max_iter {
let fx = f(x);
if fx.abs() < tol {
return x;
}
let dfx = df(x);
if dfx.abs() < 1e-50 {
break;
}
x -= fx / dfx;
}
x
}
pub fn minimize_golden_section(
f: &dyn Fn(f64) -> f64,
mut a: f64,
mut b: f64,
tol: f64,
max_iter: usize,
) -> f64 {
let phi = (5.0_f64.sqrt() - 1.0) / 2.0;
let mut c = b - phi * (b - a);
let mut d = a + phi * (b - a);
for _ in 0..max_iter {
if (b - a).abs() < tol {
break;
}
if f(c) < f(d) {
b = d;
} else {
a = c;
}
c = b - phi * (b - a);
d = a + phi * (b - a);
}
0.5 * (a + b)
}
pub fn cubic_spline_eval(x_data: &[f64], y_data: &[f64], x: f64) -> f64 {
let n = x_data.len();
if n < 2 {
return if n == 1 { y_data[0] } else { 0.0 };
}
let mut lo = 0;
let mut hi = n - 1;
while hi - lo > 1 {
let mid = (lo + hi) / 2;
if x_data[mid] > x {
hi = mid;
} else {
lo = mid;
}
}
let h = x_data[hi] - x_data[lo];
if h.abs() < 1e-50 {
return y_data[lo];
}
let a_coeff = (x_data[hi] - x) / h;
let b_coeff = (x - x_data[lo]) / h;
a_coeff * y_data[lo] + b_coeff * y_data[hi]
}
pub fn trapezoid_integrate(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n_steps: usize) -> f64 {
let h = (b - a) / n_steps as f64;
let mut sum = 0.5 * (f(a) + f(b));
for i in 1..n_steps {
let x = a + i as f64 * h;
sum += f(x);
}
sum * h
}
pub fn simpson_integrate(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n_steps: usize) -> f64 {
let n = if n_steps.is_multiple_of(2) {
n_steps
} else {
n_steps + 1
};
let h = (b - a) / n as f64;
let mut sum = f(a) + f(b);
for i in 1..n {
let x = a + i as f64 * h;
if i % 2 == 0 {
sum += 2.0 * f(x);
} else {
sum += 4.0 * f(x);
}
}
sum * h / 3.0
}
pub fn log_integrate(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n_steps: usize) -> f64 {
let ln_a = a.ln();
let ln_b = b.ln();
let d_ln = (ln_b - ln_a) / n_steps as f64;
let mut sum = 0.0;
for i in 0..n_steps {
let ln_x = ln_a + (i as f64 + 0.5) * d_ln;
let x = ln_x.exp();
sum += f(x) * x;
}
sum * d_ln
}