1pub fn shooting_method(
2 f: impl Fn(f64, &[f64]) -> Vec<f64>,
3 t_span: (f64, f64),
4 ya: f64,
5 yb: f64,
6 dt: f64,
7 tol: f64,
8 max_iter: usize,
9) -> Vec<(f64, f64)> {
10 let mut s0 = (yb - ya) / (t_span.1 - t_span.0);
11 let mut s1 = s0 + 1.0;
12
13 let integrate = |s: f64| -> f64 {
14 let mut t = t_span.0;
15 let mut y = vec![ya, s];
16 while t < t_span.1 - 1e-12 {
17 let h = dt.min(t_span.1 - t);
18 let k1 = f(t, &y);
19 let y2: Vec<f64> = y
20 .iter()
21 .zip(&k1)
22 .map(|(yi, ki)| yi + 0.5 * h * ki)
23 .collect();
24 let k2 = f(t + 0.5 * h, &y2);
25 let y3: Vec<f64> = y
26 .iter()
27 .zip(&k2)
28 .map(|(yi, ki)| yi + 0.5 * h * ki)
29 .collect();
30 let k3 = f(t + 0.5 * h, &y3);
31 let y4: Vec<f64> = y.iter().zip(&k3).map(|(yi, ki)| yi + h * ki).collect();
32 let k4 = f(t + h, &y4);
33 for (i, yi) in y.iter_mut().enumerate() {
34 *yi += h / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
35 }
36 t += h;
37 }
38 y[0]
39 };
40
41 let mut f0 = integrate(s0) - yb;
42 for _ in 0..max_iter {
43 let f1 = integrate(s1) - yb;
44 if f1.abs() < tol {
45 let mut t = t_span.0;
46 let mut y = vec![ya, s1];
47 let mut result = vec![(t, y[0])];
48 while t < t_span.1 - 1e-12 {
49 let h = dt.min(t_span.1 - t);
50 let k1 = f(t, &y);
51 let y2: Vec<f64> = y
52 .iter()
53 .zip(&k1)
54 .map(|(yi, ki)| yi + 0.5 * h * ki)
55 .collect();
56 let k2 = f(t + 0.5 * h, &y2);
57 let y3: Vec<f64> = y
58 .iter()
59 .zip(&k2)
60 .map(|(yi, ki)| yi + 0.5 * h * ki)
61 .collect();
62 let k3 = f(t + 0.5 * h, &y3);
63 let y4: Vec<f64> = y.iter().zip(&k3).map(|(yi, ki)| yi + h * ki).collect();
64 let k4 = f(t + h, &y4);
65 for (i, yi) in y.iter_mut().enumerate() {
66 *yi += h / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
67 }
68 t += h;
69 result.push((t, y[0]));
70 }
71 return result;
72 }
73 if (f1 - f0).abs() < 1e-30 {
74 break;
75 }
76 let s2 = s1 - f1 * (s1 - s0) / (f1 - f0);
77 s0 = s1;
78 f0 = f1;
79 s1 = s2;
80 }
81 vec![]
82}
83
84pub fn finite_difference_bvp(
85 p: impl Fn(f64) -> f64,
86 q: impl Fn(f64) -> f64,
87 r: impl Fn(f64) -> f64,
88 x_range: (f64, f64),
89 ya: f64,
90 yb: f64,
91 n: usize,
92) -> Vec<(f64, f64)> {
93 let h = (x_range.1 - x_range.0) / (n + 1) as f64;
94 let m = n;
95 let mut a_mat = vec![vec![0.0; m]; m];
96 let mut b_vec = vec![0.0; m];
97
98 for i in 0..m {
99 let x = x_range.0 + (i + 1) as f64 * h;
100 let pi = p(x);
101 let qi = q(x);
102 let ri = r(x);
103
104 a_mat[i][i] = -2.0 + h * h * qi;
105 if i > 0 {
106 a_mat[i][i - 1] = 1.0 - 0.5 * h * pi;
107 }
108 if i < m - 1 {
109 a_mat[i][i + 1] = 1.0 + 0.5 * h * pi;
110 }
111
112 b_vec[i] = h * h * ri;
113 if i == 0 {
114 b_vec[i] -= (1.0 - 0.5 * h * pi) * ya;
115 }
116 if i == m - 1 {
117 b_vec[i] -= (1.0 + 0.5 * h * pi) * yb;
118 }
119 }
120
121 let y_internal = tridiag_solve(&a_mat, &b_vec);
122
123 let mut result = vec![(x_range.0, ya)];
124 for (i, &yi) in y_internal.iter().enumerate() {
125 result.push((x_range.0 + (i + 1) as f64 * h, yi));
126 }
127 result.push((x_range.1, yb));
128 result
129}
130
131pub fn collocation_bvp(
132 f: impl Fn(f64, f64, f64) -> f64,
133 x_range: (f64, f64),
134 ya: f64,
135 yb: f64,
136 n: usize,
137 max_iter: usize,
138 tol: f64,
139) -> Vec<(f64, f64)> {
140 let h = (x_range.1 - x_range.0) / (n + 1) as f64;
141 let mut y = vec![0.0; n];
142 for (i, yi) in y.iter_mut().enumerate() {
143 let t = (i + 1) as f64 / (n + 1) as f64;
144 *yi = ya * (1.0 - t) + yb * t;
145 }
146 for _ in 0..max_iter {
147 let mut residual = vec![0.0; n];
148 for i in 0..n {
149 let x = x_range.0 + (i + 1) as f64 * h;
150 let y_prev = if i == 0 { ya } else { y[i - 1] };
151 let y_next = if i == n - 1 { yb } else { y[i + 1] };
152 let ypp = (y_prev - 2.0 * y[i] + y_next) / (h * h);
153 let yp = (y_next - y_prev) / (2.0 * h);
154 residual[i] = ypp - f(x, y[i], yp);
155 }
156 let max_res: f64 = residual.iter().map(|r| r.abs()).fold(0.0, f64::max);
157 if max_res < tol {
158 break;
159 }
160 for (yi, &ri) in y.iter_mut().zip(residual.iter()) {
161 *yi += 0.5 * h * h * ri;
162 }
163 }
164 let mut result = vec![(x_range.0, ya)];
165 for (i, &yi) in y.iter().enumerate() {
166 result.push((x_range.0 + (i + 1) as f64 * h, yi));
167 }
168 result.push((x_range.1, yb));
169 result
170}
171
172pub fn multiple_shooting(
173 f: impl Fn(f64, &[f64]) -> Vec<f64>,
174 t_span: (f64, f64),
175 ya: f64,
176 yb: f64,
177 segments: usize,
178 dt: f64,
179 tol: f64,
180 max_iter: usize,
181) -> Vec<(f64, f64)> {
182 let seg_len = (t_span.1 - t_span.0) / segments as f64;
183 let mut slopes = vec![(yb - ya) / (t_span.1 - t_span.0); segments];
184 for _ in 0..max_iter {
185 let mut all_good = true;
186 for seg in 0..segments {
187 let t0 = t_span.0 + seg as f64 * seg_len;
188 let t1 = t0 + seg_len;
189 let y0_seg = if seg == 0 {
190 ya
191 } else {
192 ya + slopes[..seg].iter().sum::<f64>() * seg_len
193 };
194 let mut t = t0;
195 let mut y = vec![y0_seg, slopes[seg]];
196 while t < t1 - 1e-12 {
197 let h = dt.min(t1 - t);
198 let k1 = f(t, &y);
199 let y2: Vec<f64> = y
200 .iter()
201 .zip(&k1)
202 .map(|(yi, ki)| yi + 0.5 * h * ki)
203 .collect();
204 let k2 = f(t + 0.5 * h, &y2);
205 let y3: Vec<f64> = y
206 .iter()
207 .zip(&k2)
208 .map(|(yi, ki)| yi + 0.5 * h * ki)
209 .collect();
210 let k3 = f(t + 0.5 * h, &y3);
211 let y4: Vec<f64> = y.iter().zip(&k3).map(|(yi, ki)| yi + h * ki).collect();
212 let k4 = f(t + h, &y4);
213 for (i, yi) in y.iter_mut().enumerate() {
214 *yi += h / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
215 }
216 t += h;
217 }
218 let target = if seg == segments - 1 {
219 yb
220 } else {
221 ya + slopes[..=seg].iter().sum::<f64>() * seg_len
222 };
223 let err = y[0] - target;
224 if err.abs() > tol {
225 all_good = false;
226 }
227 slopes[seg] -= err / seg_len;
228 }
229 if all_good {
230 break;
231 }
232 }
233 shooting_method(&f, t_span, ya, yb, dt, tol, max_iter)
234}
235
236pub fn sturm_liouville_eigenvalues(
237 p: impl Fn(f64) -> f64,
238 q: impl Fn(f64) -> f64,
239 w: impl Fn(f64) -> f64,
240 x_range: (f64, f64),
241 n: usize,
242 num_eigenvalues: usize,
243) -> Vec<f64> {
244 let h = (x_range.1 - x_range.0) / (n + 1) as f64;
245 let m = n;
246 let mut a_mat = vec![vec![0.0; m]; m];
247 let mut b_mat = vec![vec![0.0; m]; m];
248 for i in 0..m {
249 let x = x_range.0 + (i + 1) as f64 * h;
250 let pi = p(x);
251 let qi = q(x);
252 let wi = w(x);
253 a_mat[i][i] = 2.0 * pi / (h * h) + qi;
254 if i > 0 {
255 a_mat[i][i - 1] = -pi / (h * h);
256 }
257 if i < m - 1 {
258 a_mat[i][i + 1] = -pi / (h * h);
259 }
260 b_mat[i][i] = wi;
261 }
262 let mut eigenvalues = Vec::new();
263 let mut lambda = 0.1;
264 for _ in 0..num_eigenvalues {
265 for _ in 0..100 {
266 let mut mat = vec![vec![0.0; m]; m];
267 for (i, mat_i) in mat.iter_mut().enumerate() {
268 for (j, mat_ij) in mat_i.iter_mut().enumerate() {
269 *mat_ij = a_mat[i][j] - lambda * b_mat[i][j];
270 }
271 }
272 let det = matrix_det_small(&mat);
273 let eps = 1e-6;
274 let mut mat2 = mat.clone();
275 for (i, m2i) in mat2.iter_mut().enumerate() {
276 m2i[i] -= eps;
277 }
278 let det2 = matrix_det_small(&mat2);
279 let ddet = (det2 - det) / (-eps);
280 if ddet.abs() < 1e-30 {
281 break;
282 }
283 let correction = det / ddet;
284 lambda -= correction;
285 if correction.abs() < 1e-10 {
286 break;
287 }
288 }
289 eigenvalues.push(lambda);
290 lambda += 1.0;
291 }
292 eigenvalues
293}
294
295fn matrix_det_small(m: &[Vec<f64>]) -> f64 {
296 let n = m.len();
297 let mut lu = m.to_vec();
298 let mut det = 1.0;
299 for j in 0..n {
300 for i in j + 1..n {
301 if lu[j][j].abs() < 1e-30 {
302 return 0.0;
303 }
304 let factor = lu[i][j] / lu[j][j];
305 let (top, bot) = lu.split_at_mut(i);
306 for (d, &s) in bot[0][j..n].iter_mut().zip(&top[j][j..n]) {
307 *d -= factor * s;
308 }
309 }
310 det *= lu[j][j];
311 }
312 det
313}
314
315pub fn relaxation_bvp(
316 f: impl Fn(f64, f64) -> f64,
317 x_range: (f64, f64),
318 ya: f64,
319 yb: f64,
320 n: usize,
321 omega: f64,
322 max_iter: usize,
323 tol: f64,
324) -> Vec<(f64, f64)> {
325 let h = (x_range.1 - x_range.0) / (n + 1) as f64;
326 let mut y = vec![0.0; n + 2];
327 y[0] = ya;
328 y[n + 1] = yb;
329 for (i, yi) in y.iter_mut().enumerate().skip(1).take(n) {
330 let t = i as f64 / (n + 1) as f64;
331 *yi = ya * (1.0 - t) + yb * t;
332 }
333 for _ in 0..max_iter {
334 let mut max_change = 0.0_f64;
335 for i in 1..=n {
336 let x = x_range.0 + i as f64 * h;
337 let y_new = 0.5 * (y[i - 1] + y[i + 1] - h * h * f(x, y[i]));
338 let change = (y_new - y[i]).abs();
339 max_change = max_change.max(change);
340 y[i] = (1.0 - omega) * y[i] + omega * y_new;
341 }
342 if max_change < tol {
343 break;
344 }
345 }
346 (0..=n + 1)
347 .map(|i| (x_range.0 + i as f64 * h, y[i]))
348 .collect()
349}
350
351fn tridiag_solve(a: &[Vec<f64>], b: &[f64]) -> Vec<f64> {
352 let n = b.len();
353 let mut c = vec![0.0; n];
354 let mut d = b.to_vec();
355
356 c[0] = a[0][1] / a[0][0];
357 d[0] /= a[0][0];
358
359 for i in 1..n {
360 let sub = if i > 0 { a[i][i - 1] } else { 0.0 };
361 let sup = if i < n - 1 { a[i][i + 1] } else { 0.0 };
362 let m = a[i][i] - sub * c[i - 1];
363 if i < n - 1 {
364 c[i] = sup / m;
365 }
366 d[i] = (d[i] - sub * d[i - 1]) / m;
367 }
368
369 let mut x = vec![0.0; n];
370 x[n - 1] = d[n - 1];
371 for i in (0..n - 1).rev() {
372 x[i] = d[i] - c[i] * x[i + 1];
373 }
374 x
375}