Open MucTepDayH16 opened 1 year ago
Hi. Unfortunately I won't have time to work on this in the near future. Feel free to implement it and submit a PR.
You can add support for complex vectors in Rk4
by just removing some unused generic constraints and adding T: From<f64>
for multiplying the vectors by f64
s (see this git diff).
This seems to work fine for Rk4
, but the removed constraints are used by the other solvers. I'm not sure if they can be changed in the same way, and if not, they will have different constraints than Rk4
.
Support for Complex ode systems could be implemented by wrapping the existing real solvers, similar to the approach used in scipy with the complex_ode class
The mathematical background is that the complex field $\mathbb C$ can be seen as a linear combination over a two dimensional real field $\mathbb R^2$. The basis is formed by 1 and the imaginary unit (i as solution to the polynomial $X^ 2+1=0$). An n dimensional complex state vector can thus be seen as a 2n dimensional real state vector by splitting the real and imaginary parts.
Instead of directly solving the complex ode system, we would map the n dimensional complex system to a 2n dimensional real system. Then, we can apply all the existing (and future) real ode solvers to the real state, and map the results back to the complex state.
The user defines the system using complex numbers, i.e., with nalgebra::Complex An implementation could wrap every call to the ˋsystemˋ method of the user-defined complex ˋSystemˋ with the mapping of complex to real state vectors. Internally, the ode-solvers then only consider real-valued vectors. @srenevey what do you think?
The exponential equation has the form $$\frac{dz}{dt}=kz $$ with $k\in\mathbb C$ being a complex constant and $z:\mathbb R \to\mathbb C$ is a complex function of a real variable $t$ (e.g. time).
The solution to this equation is $$z(t)=z(0)e^{kt}$$ wherein $z(0)$ is the initial value at the begin of the time interval. Using for instance $k=i$ and $z(0)=1$, the solution is $z(t)=e^{it}=\cos t + i \sin t$ (Eulers Formula). See MIT sup 6 PDF for additional information.
Instead of solving the complex system for $z\in\mathbb C$, we solve for the real coefficients in $z=a+bi$ with $(a,b)\in\mathbb R^2$. In our example this gives $a=cos(t), b=sin(t)$.
#[cfg(test)]
mod tests {
use crate::rk4::Rk4;
use crate::{OVector, System, Vector2};
use nalgebra::{allocator::Allocator, DefaultAllocator, Dim};
struct Test1 {}
impl<D: Dim> System<f64, OVector<f64, D>> for Test1
where
DefaultAllocator: Allocator<f64, D>,
{
fn system(&self, x: f64, y: &OVector<f64, D>, dy: &mut OVector<f64, D>) {
// dy/dt = iy with y = a + bi => dy/dt = ai - b = -b + ai
dy[0] = -y[1]; // -b
dy[1] = y[0]; // a
}
}
#[test]
fn test_complex_test1_svector() {
let system = Test1 {};
let t_end: f64 = 0.2;
let real_end_state = Vector2::new(t_end.cos(), t_end.sin());
let mut stepper = Rk4::new(system, 0., Vector2::new(1., 0.), t_end, 0.1);
let _ = stepper.integrate();
let x_out = stepper.x_out();
let y_out = stepper.y_out();
assert!(
(*x_out.last().unwrap() - t_end).abs() < 1.0E-8,
"x_out must end with x_end"
);
assert!(
y_out
.last()
.unwrap()
.relative_eq(&real_end_state, 1.0E-5, 1.0E-5),
"The last state must be the solution at x_end!"
);
}
}
I've implemented
System
trait for my System, but I want to useDVector<Complex<f64>>
as myy
value.As i can see,
DVector
support in this crate is WIP, but nothing is said aboutComplex
types. Would you manage to add support for it?