1pub fn lotka_volterra(
2 alpha: f64,
3 beta: f64,
4 delta: f64,
5 gamma: f64,
6) -> impl Fn(f64, &[f64]) -> Vec<f64> {
7 move |_t: f64, y: &[f64]| {
8 vec![
9 alpha * y[0] - beta * y[0] * y[1],
10 delta * y[0] * y[1] - gamma * y[1],
11 ]
12 }
13}
14
15pub fn lorenz(sigma: f64, rho: f64, beta: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
16 move |_t: f64, y: &[f64]| {
17 vec![
18 sigma * (y[1] - y[0]),
19 y[0] * (rho - y[2]) - y[1],
20 y[0] * y[1] - beta * y[2],
21 ]
22 }
23}
24
25pub fn van_der_pol(mu: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
26 move |_t: f64, y: &[f64]| vec![y[1], mu * (1.0 - y[0] * y[0]) * y[1] - y[0]]
27}
28
29pub fn harmonic_oscillator(omega: f64, damping: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
30 move |_t: f64, y: &[f64]| vec![y[1], -2.0 * damping * omega * y[1] - omega * omega * y[0]]
31}
32
33pub fn double_pendulum(
34 l1: f64,
35 l2: f64,
36 m1: f64,
37 m2: f64,
38 g: f64,
39) -> impl Fn(f64, &[f64]) -> Vec<f64> {
40 move |_t: f64, y: &[f64]| {
41 let (th1, w1, th2, w2) = (y[0], y[1], y[2], y[3]);
42 let delta = th1 - th2;
43 let den1 = (m1 + m2) * l1 - m2 * l1 * delta.cos() * delta.cos();
44 let den2 = (l2 / l1) * den1;
45
46 let dw1 = (m2 * l1 * w1 * w1 * delta.sin() * delta.cos()
47 + m2 * g * th2.sin() * delta.cos()
48 + m2 * l2 * w2 * w2 * delta.sin()
49 - (m1 + m2) * g * th1.sin())
50 / den1;
51
52 let dw2 = (-m2 * l2 * w2 * w2 * delta.sin() * delta.cos()
53 + (m1 + m2) * g * th1.sin() * delta.cos()
54 - (m1 + m2) * l1 * w1 * w1 * delta.sin()
55 - (m1 + m2) * g * th2.sin())
56 / den2;
57
58 vec![w1, dw1, w2, dw2]
59 }
60}
61
62pub fn sir_model(beta: f64, gamma: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
63 move |_t: f64, y: &[f64]| {
64 let (s, i, _) = (y[0], y[1], y[2]);
65 vec![-beta * s * i, beta * s * i - gamma * i, gamma * i]
66 }
67}
68
69pub fn rossler(a: f64, b: f64, c: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
70 move |_t: f64, y: &[f64]| vec![-y[1] - y[2], y[0] + a * y[1], b + y[2] * (y[0] - c)]
71}
72
73pub fn three_body_planar(m: [f64; 3], g: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
74 move |_t: f64, y: &[f64]| {
75 let mut dy = vec![0.0; 12];
76 for i in 0..3 {
77 dy[2 * i] = y[6 + 2 * i];
78 dy[2 * i + 1] = y[6 + 2 * i + 1];
79 }
80 for i in 0..3 {
81 for j in 0..3 {
82 if i == j {
83 continue;
84 }
85 let dx = y[2 * j] - y[2 * i];
86 let dy_val = y[2 * j + 1] - y[2 * i + 1];
87 let r = (dx * dx + dy_val * dy_val).sqrt();
88 if r < 1e-30 {
89 continue;
90 }
91 let r3 = r * r * r;
92 dy[6 + 2 * i] += g * m[j] * dx / r3;
93 dy[6 + 2 * i + 1] += g * m[j] * dy_val / r3;
94 }
95 }
96 dy
97 }
98}
99
100pub fn brusselator(a: f64, b: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
101 move |_t: f64, y: &[f64]| {
102 vec![
103 a - (b + 1.0) * y[0] + y[0] * y[0] * y[1],
104 b * y[0] - y[0] * y[0] * y[1],
105 ]
106 }
107}
108
109pub fn oregonator(epsilon: f64, q: f64, f_param: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
110 move |_t: f64, y: &[f64]| {
111 vec![
112 (q * y[1] - y[0] * y[1] + y[0] * (1.0 - y[0])) / epsilon,
113 -q * y[1] - y[0] * y[1] + f_param * y[2],
114 y[0] - y[2],
115 ]
116 }
117}
118
119pub fn fitzhugh_nagumo(a: f64, b: f64, tau: f64, i_ext: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
120 move |_t: f64, y: &[f64]| {
121 vec![
122 y[0] - y[0] * y[0] * y[0] / 3.0 - y[1] + i_ext,
123 (y[0] + a - b * y[1]) / tau,
124 ]
125 }
126}
127
128pub fn duffing(
129 alpha: f64,
130 beta: f64,
131 gamma: f64,
132 delta: f64,
133 omega: f64,
134) -> impl Fn(f64, &[f64]) -> Vec<f64> {
135 move |t: f64, y: &[f64]| {
136 vec![
137 y[1],
138 -delta * y[1] - alpha * y[0] - beta * y[0].powi(3) + gamma * (omega * t).cos(),
139 ]
140 }
141}
142
143pub fn chen_system(a: f64, b: f64, c: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
144 move |_t: f64, y: &[f64]| {
145 vec![
146 a * (y[1] - y[0]),
147 (c - a) * y[0] - y[0] * y[2] + c * y[1],
148 y[0] * y[1] - b * y[2],
149 ]
150 }
151}
152
153pub fn chua_circuit(alpha: f64, beta: f64, m0: f64, m1: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
154 move |_t: f64, y: &[f64]| {
155 let h = m1 * y[0] + 0.5 * (m0 - m1) * ((y[0] + 1.0).abs() - (y[0] - 1.0).abs());
156 vec![alpha * (y[1] - y[0] - h), y[0] - y[1] + y[2], -beta * y[1]]
157 }
158}
159
160pub fn predator_prey_holling(
161 r: f64,
162 k: f64,
163 a: f64,
164 b: f64,
165 c: f64,
166 d: f64,
167) -> impl Fn(f64, &[f64]) -> Vec<f64> {
168 move |_t: f64, y: &[f64]| {
169 let prey = y[0];
170 let pred = y[1];
171 let functional_response = a * prey / (1.0 + b * prey);
172 vec![
173 r * prey * (1.0 - prey / k) - functional_response * pred,
174 c * functional_response * pred - d * pred,
175 ]
176 }
177}
178
179pub fn stiff_robertson(k1: f64, k2: f64, k3: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
180 move |_t: f64, y: &[f64]| {
181 vec![
182 -k1 * y[0] + k3 * y[1] * y[2],
183 k1 * y[0] - k2 * y[1] * y[1] - k3 * y[1] * y[2],
184 k2 * y[1] * y[1],
185 ]
186 }
187}
188
189pub fn rigid_body(i1: f64, i2: f64, i3: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
190 move |_t: f64, y: &[f64]| {
191 vec![
192 (i2 - i3) / i1 * y[1] * y[2],
193 (i3 - i1) / i2 * y[0] * y[2],
194 (i1 - i2) / i3 * y[0] * y[1],
195 ]
196 }
197}
198
199pub fn seir_model(beta: f64, sigma: f64, gamma: f64) -> impl Fn(f64, &[f64]) -> Vec<f64> {
200 move |_t: f64, y: &[f64]| {
201 let (s, e, i, _) = (y[0], y[1], y[2], y[3]);
202 vec![
203 -beta * s * i,
204 beta * s * i - sigma * e,
205 sigma * e - gamma * i,
206 gamma * i,
207 ]
208 }
209}