1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

use crate::algebra::linear::{Vector};
use crate::algebra::abstr::Real;
use super::explicit_method::{ExplicitFixedStepSizeMethod};
use super::ExplicitODE;
use crate::analysis::ode::fixed_stepper::ExplicitFixedStepper;

/// Solves an ordinary differential equation using the 3th order Runge-Kutta algorithm.
///
///<a href="https://en.wikipedia.org/wiki/Runge-Kutta_methods">https://en.wikipedia
/// .org/wiki/Rung-Kutta_methods</a>
pub struct Kutta3<T>
{
    stepper: ExplicitFixedStepper<T>
}

impl<T> Kutta3<T>
    where T: Real
{
    /// Creates a Kutta3 instance with step size 'step_size'
    ///
    pub fn new(step_size: T) -> Kutta3<T>
    {
        return Kutta3
        {
            stepper: ExplicitFixedStepper::new(step_size)
        }
    }

    pub fn solve<F>(self: &Self, prob: &F) -> Result<(Vec<T>, Vec<Vector<T>>), ()>
        where F: ExplicitODE<T>
    {
        return self.stepper.solve(prob, self);
    }

    pub fn get_step_size(self: &Self) -> &T
    {
        return self.stepper.get_step_size();
    }

    pub fn set_step_size(self: &mut Self, step_size: T)
    {
        self.stepper.set_step_size(step_size)
    }
}

impl<T> ExplicitFixedStepSizeMethod<T> for Kutta3<T>
    where T: Real
{
    fn do_step<F>(self: &Self, prob: &F, t_n: &T, x_n: &Vector<T>, h: &T) -> Vector<T>
        where F: ExplicitODE<T>
    {
        // k1 = f(t_n, x_n)
        let k1: Vector<T> = prob.func(&t_n, &x_n);

        // k2 = func(t_n + h / 2, x_n + h / 2 k1)
        let k2: Vector<T> = prob.func(&(*t_n + (*h / T::from_f64(2.0).unwrap())), &(x_n + &((&k1 * h) / T::from_f64(2.0).unwrap())));

        // k3 = func(t_n + h, x_n + - h k1 + 2h *k2)
        let k3: Vector<T> = prob.func(&(*t_n + *h), &(x_n + &(&(&(&k2 * &T::from_f64(2.0).unwrap()) - &k1) * h))) ;

        // x[n+1] = xn + h*(k1 + 4*k2 + k3)/6
        return x_n + &((k1 + k2 * T::from_f64(4.0).unwrap() + k3) * *h / T::from_f64(6.0).unwrap());
    }

    /// Kuttas method is a 3rd order method
    fn order(self: &Self) -> u8
    {
        return 3;
    }
}