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

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 Heun's method.
///
/// <a href="https://en.wikipedia.org/wiki/Heun's_method">https://en.wikipedia.org/wiki/Heun's_method</a>
pub struct Heun<T>
{
    stepper: ExplicitFixedStepper<T>
}


impl<T> Heun<T>
    where T: Real
{
    /// Creates a Heun instance with step size 'step_size'
    ///
    pub fn new(step_size: T) -> Heun<T>
    {
        return Heun
        {
            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 Heun<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>
     {
            let k_1: Vector<T> = prob.func(t_n, x_n);
            let k_2: Vector<T> = prob.func(&(*t_n + *h), &(x_n + &(k_1.clone() * *h)));

            return x_n.clone() + (&k_1 + &k_2) * *h / T::from_f64(2.0);
    }

    ///
    fn order(self: &Self) -> u8
    {
        return 2;
    }
}