scirs2_integrate/dae/utils/
mod.rs

1//! Utility functions for DAE solvers
2//!
3//! This module provides utility functions for use in DAE solvers.
4
5use crate::IntegrateFloat;
6use scirs2_core::ndarray::{Array2, ArrayView2};
7// use scirs2_core::numeric::{Float, FromPrimitive};
8
9// Linear solvers
10pub mod linear_solvers;
11
12// Re-export useful utilities
13pub use linear_solvers::{matrix_norm, solve_linear_system, solve_lu, vector_norm};
14
15/// Compute the constraint Jacobian for a constraint function
16#[allow(dead_code)]
17pub fn compute_constraint_jacobian<F: IntegrateFloat>(
18    g: &impl Fn(F, &[F], &[F]) -> Vec<F>,
19    t: F,
20    x: &[F],
21    y: &[F],
22) -> Array2<F> {
23    let n = x.len();
24    let m = y.len();
25    let epsilon = F::from_f64(1e-8).unwrap();
26
27    // Compute g at the current point
28    let g0 = g(t, x, y);
29    let ng = g0.len();
30
31    let mut jacobian = Array2::zeros((ng, n + m));
32
33    // Compute partial derivatives with respect to x
34    for i in 0..n {
35        let mut x_perturbed = x.to_vec();
36        x_perturbed[i] += epsilon;
37        let g_perturbed = g(t, &x_perturbed, y);
38
39        for j in 0..ng {
40            jacobian[(j, i)] = (g_perturbed[j] - g0[j]) / epsilon;
41        }
42    }
43
44    // Compute partial derivatives with respect to y
45    for i in 0..m {
46        let mut y_perturbed = y.to_vec();
47        y_perturbed[i] += epsilon;
48        let g_perturbed = g(t, x, &y_perturbed);
49
50        for j in 0..ng {
51            jacobian[(j, n + i)] = (g_perturbed[j] - g0[j]) / epsilon;
52        }
53    }
54
55    jacobian
56}
57
58/// Check if a matrix is singular
59#[allow(dead_code)]
60pub fn is_singular_matrix<F: IntegrateFloat>(matrix: ArrayView2<F>) -> bool {
61    let n = matrix.nrows();
62    if n != matrix.ncols() {
63        return true; // Non-square matrices are considered singular
64    }
65
66    // Use a simple determinant check for small matrices
67    if n == 1 {
68        return matrix[(0, 0)].abs() < F::from_f64(1e-10).unwrap();
69    }
70
71    // For larger matrices, check condition number or use LU decomposition
72    // For now, use a simple diagonal dominance check as a heuristic
73    let epsilon = F::from_f64(1e-10).unwrap();
74
75    for i in 0..n {
76        let diagonal = matrix[(i, i)].abs();
77        let mut off_diagonal_sum = F::zero();
78
79        for j in 0..n {
80            if i != j {
81                off_diagonal_sum += matrix[(i, j)].abs();
82            }
83        }
84
85        if diagonal < epsilon || diagonal < off_diagonal_sum * F::from_f64(0.1).unwrap() {
86            return true;
87        }
88    }
89
90    false
91}
scirs2_integrate/dae/utils/mod.rs

scirs2_integrate/dae/utils/
mod.rs
