Implement system solvers for complex vectors

[d-drozhzhin]

I've implemented System trait for my System, but I want to use DVector<Complex<f64>> as my y value.

As i can see, DVector support in this crate is WIP, but nothing is said about Complex types. Would you manage to add support for it?

[srenevey]

Hi. Unfortunately I won't have time to work on this in the near future. Feel free to implement it and submit a PR.

[jonah-h]

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 f64s (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.

[tomtuamnuq]

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?

Example

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)$.

Unit Test (draft with real system)

#[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!"
        );
    }
}