harmonic_oscillator/
harmonic_oscillator.rs

1
2extern crate maths_traits;
3extern crate numerical_integration;
4
5use numerical_integration::{VelocityVerlet, VelIntegrator};
6
7
8fn main() {
9
10    use maths_traits::algebra::*;
11    // use maths_traits::analysis::*;
12    // use maths_traits::analysis::Exponential;
13
14    //
15    //A simple 2D vector impl
16    //
17
18    #[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)]
19    pub struct Vec2<T> {
20        pub x: T,
21        pub y: T,
22    }
23
24    // impl<T> Vec2<T> {
25    //     #[inline]
26    //     pub fn new(x: T, y: T) -> Self {
27    //         Vec2 { x: x, y: y }
28    //     }
29    // }
30
31    impl<T: Add<Output=T>> Add for Vec2<T> {
32        type Output = Self;
33        #[inline]
34        fn add(self, rhs: Self) -> Self {
35            Vec2 {
36                x: self.x + rhs.x,
37                y: self.y + rhs.y,
38            }
39        }
40    }
41
42    impl<T: AddAssign> AddAssign for Vec2<T> {
43        #[inline]
44        fn add_assign(&mut self, rhs: Self) {
45            self.x += rhs.x;
46            self.y += rhs.y;
47        }
48    }
49
50    impl<T: Sub<Output = T>> Sub for Vec2<T> {
51        type Output = Self;
52        #[inline]
53        fn sub(self, rhs: Self) -> Self {
54            Vec2 {
55                x: self.x - rhs.x,
56                y: self.y - rhs.y,
57            }
58        }
59    }
60
61    impl<T: SubAssign> SubAssign for Vec2<T> {
62        #[inline]
63        fn sub_assign(&mut self, rhs: Self) {
64            self.x -= rhs.x;
65            self.y -= rhs.y;
66        }
67    }
68
69    impl<T: Neg<Output = T>> Neg for Vec2<T> {
70        type Output = Self;
71        #[inline]
72        fn neg(self) -> Self {
73            Vec2 {
74                x: -self.x,
75                y: -self.y,
76            }
77        }
78    }
79
80    impl<T: Zero> Zero for Vec2<T> {
81        #[inline]
82        fn zero() -> Self {
83            Vec2 {
84                x: T::zero(),
85                y: T::zero()
86            }
87        }
88        #[inline]
89        fn is_zero(&self) -> bool {
90            self.x.is_zero() && self.y.is_zero()
91        }
92    }
93
94    impl<K: Clone, T: Mul<K, Output = T>> Mul<K> for Vec2<T> {
95        type Output = Self;
96        #[inline]
97        fn mul(self, rhs: K) -> Self {
98            Vec2 {
99                x: self.x * rhs.clone(),
100                y: self.y * rhs,
101            }
102        }
103    }
104
105    impl<K: Clone, T: Clone + MulAssign<K>> MulAssign<K> for Vec2<T> {
106        #[inline]
107        fn mul_assign(&mut self, rhs: K) {
108            self.x *= rhs.clone();
109            self.y *= rhs;
110        }
111    }
112
113    impl<K: Clone, T: Div<K, Output = T>> Div<K> for Vec2<T> {
114        type Output = Self;
115        #[inline]
116        fn div(self, rhs: K) -> Self {
117            Vec2 {
118                x: self.x / rhs.clone(),
119                y: self.y / rhs,
120            }
121        }
122    }
123
124    impl<K: Clone, T: Clone + DivAssign<K>> DivAssign<K> for Vec2<T> {
125        #[inline]
126        fn div_assign(&mut self, rhs: K) {
127            self.x /= rhs.clone();
128            self.y /= rhs;
129        }
130    }
131
132    impl<T> AddAssociative for Vec2<T> {}
133    impl<T> MulAssociative for Vec2<T> {}
134    impl<T> AddCommutative for Vec2<T> {}
135    impl<T> MulCommutative for Vec2<T> {}
136    impl<T> Distributive<T> for Vec2<T> {}
137
138    //
139    //A simple harmonic oscillator simulation
140    //
141    //we use a vec2 for the state vector, where the x coord is the position
142    //and the y coord is the velocity
143    //
144
145    //the x-coord of the derivative is y since it is the velocity
146    //and the y coord of the derivative is -x as dv/dt = -x
147    fn f(_t: f64, y: Vec2<f64>) -> Vec2<f64> { Vec2 { x:y.y, y:-y.x } }
148
149    //shifts the velocity into x and sets y to be 0
150    fn v(_t: f64, y: Vec2<f64>) -> Vec2<f64> { Vec2 { x:y.y, y:0.0 } }
151
152    //the time-step
153    let dt = 0.1;
154
155    //init
156    let mut t = 0.0;
157    let y0 = Vec2{x:1.0, y:0.0};
158    let mut y1 = VelocityVerlet.init_with_vel(y0, dt, &v, &f);
159
160    // let mut y2 = RK4.init(y0, dt, &f);
161    // let mut y3 = RK4.init(y0, dt, &f);
162
163    for _ in 0..100 {
164
165        let yi = VelocityVerlet.step_with_vel(t, y1.as_mut(), dt, &v, &f);
166        t += dt;
167
168        println!(
169            "{:>6.3} | {:>10.7} {:>10.7} | {:>10.7} {:>10.6} ",
170            t, yi.x, yi.y, t.cos(), -t.sin()
171        );
172    }
173
174
175}
harmonic_oscillator/harmonic_oscillator.rs

harmonic_oscillator/
harmonic_oscillator.rs
