mathhook-core 0.2.0

//! Homogeneous ODE solver
//!
//! Solves homogeneous ODEs of the form dy/dx = f(y/x).
//!
//! A first-order ODE is homogeneous if it can be written as dy/dx = f(y/x).
//! The substitution v = y/x transforms it into a separable ODE.

use super::ODEResult;
use crate::calculus::integrals::Integration;
use crate::core::{Expression, Symbol};
use crate::simplify::Simplify;
use crate::symbol;

/// Homogeneous ODE solver
///
/// Solves ODEs of the form dy/dx = f(y/x) using substitution v = y/x.
pub struct HomogeneousODESolver;

impl HomogeneousODESolver {
    /// Solve homogeneous ODE
    ///
    /// Transforms homogeneous ODE to separable form via substitution v = y/x.
    ///
    /// # Arguments
    ///
    /// * `f` - Right-hand side function f(y/x)
    /// * `dependent` - Dependent variable (y)
    /// * `independent` - Independent variable (x)
    ///
    /// # Returns
    ///
    /// General solution y(x) or implicit solution
    ///
    /// # Examples
    ///
    /// ```rust
    /// use mathhook_core::calculus::ode::first_order::homogeneous::HomogeneousODESolver;
    /// use mathhook_core::{symbol, expr};
    ///
    /// let x = symbol!(x);
    /// let y = symbol!(y);
    /// let v = symbol!(v);
    ///
    /// // dy/dx = (y/x)^2 (homogeneous)
    /// // f(v) = v^2 where v = y/x
    /// let f = expr!(v ^ 2);
    ///
    /// let solver = HomogeneousODESolver;
    /// let solution = solver.solve(&f, &y, &x);
    /// assert!(solution.is_ok());
    /// ```
    pub fn solve(&self, f: &Expression, _dependent: &Symbol, independent: &Symbol) -> ODEResult {
        // Original equation: dy/dx = f(y/x)
        // Substitution: v = y/x, so y = vx
        // Then: dy/dx = v + x(dv/dx)

        // Substituting into original equation:
        // v + x(dv/dx) = f(v)
        // x(dv/dx) = f(v) - v
        // dv/(f(v) - v) = dx/x

        // Create v symbol for substitution
        let v = symbol!(v);

        // Compute f(v) - v
        let f_minus_v = Expression::add(vec![
            f.clone(),
            Expression::mul(vec![Expression::integer(-1), Expression::symbol(v.clone())]),
        ])
        .simplify();

        // Separate variables: dv/(f(v) - v) = dx/x
        // Integrate both sides

        // Left side: ∫ dv/(f(v) - v)
        let lhs_integrand = Expression::pow(f_minus_v, Expression::integer(-1)).simplify();
        let lhs_integral = lhs_integrand.integrate(v, 0);

        // Right side: ∫ dx/x = ln|x|
        let rhs_integral =
            Expression::function("log", vec![Expression::symbol(independent.clone())]);

        // Solution: ∫ dv/(f(v) - v) = ln|x| + C
        let c = Expression::symbol(symbol!(C));
        let solution = Expression::add(vec![
            lhs_integral,
            Expression::mul(vec![Expression::integer(-1), rhs_integral]),
            c,
        ])
        .simplify();

        Ok(solution)
    }

    /// Detect if ODE is homogeneous
    ///
    /// Checks if dy/dx can be expressed as function of y/x only.
    ///
    /// # Arguments
    ///
    /// * `equation` - The ODE equation
    /// * `dependent` - Dependent variable (y)
    /// * `independent` - Independent variable (x)
    ///
    /// # Returns
    ///
    /// True if equation is homogeneous
    pub fn is_homogeneous(
        &self,
        _equation: &Expression,
        _dependent: &Symbol,
        _independent: &Symbol,
    ) -> bool {
        // Pattern matching for homogeneous form would go here
        // This is a placeholder for future enhancement
        false
    }

    /// Extract homogeneous function f(v) from ODE
    ///
    /// Given dy/dx = expression, extract f(v) where v = y/x
    ///
    /// # Arguments
    ///
    /// * `rhs` - Right-hand side of dy/dx = rhs
    /// * `dependent` - Dependent variable (y)
    /// * `independent` - Independent variable (x)
    ///
    /// # Returns
    ///
    /// Function f(v) if successfully extracted
    pub fn extract_homogeneous_function(
        &self,
        rhs: &Expression,
        dependent: &Symbol,
        independent: &Symbol,
    ) -> Option<Expression> {
        // Replace all occurrences of y with v*x
        // Then simplify and check if result is function of v only
        // This is a placeholder for future enhancement
        let _ = (rhs, dependent, independent);
        None
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{expr, symbol};

    #[test]
    fn test_homogeneous_simple() {
        let x = symbol!(x);
        let y = symbol!(y);

        // dy/dx = y/x (f(v) = v)
        let f = expr!(v);

        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
        let sol = solution.unwrap();
        let sol_str = sol.to_string();
        // Solution should involve logarithm
        assert!(sol_str.contains("log") || sol_str.contains("ln"));
    }

    #[test]
    fn test_homogeneous_quadratic() {
        let x = symbol!(x);
        let y = symbol!(y);

        // dy/dx = (y/x)^2 = v^2
        let f = expr!(v ^ 2);

        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
    }

    #[test]
    fn test_homogeneous_linear_in_v() {
        let x = symbol!(x);
        let y = symbol!(y);

        // dy/dx = 2(y/x) = 2v
        let f = expr!(2 * v);

        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
    }

    #[test]
    fn test_homogeneous_rational() {
        let x = symbol!(x);
        let y = symbol!(y);

        // dy/dx = (y/x + 1) / (y/x - 1) = (v + 1) / (v - 1)
        let f = Expression::mul(vec![
            Expression::add(vec![expr!(v), Expression::integer(1)]),
            Expression::pow(
                Expression::add(vec![
                    expr!(v),
                    Expression::mul(vec![Expression::integer(-1), Expression::integer(1)]),
                ]),
                Expression::integer(-1),
            ),
        ]);

        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
    }

    #[test]
    fn test_homogeneous_with_trig() {
        let x = symbol!(x);
        let y = symbol!(y);

        // dy/dx = sin(y/x) = sin(v)
        let f = expr!(sin(v));

        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
    }

    #[test]
    fn test_homogeneous_solution_structure() {
        let x = symbol!(x);
        let y = symbol!(y);

        let f = expr!(v ^ 2);
        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x).unwrap();

        // Solution should contain constant C
        let sol_str = solution.to_string();
        assert!(sol_str.contains("C"));
    }

    #[test]
    fn test_homogeneous_implicit_solution() {
        let x = symbol!(x);
        let y = symbol!(y);

        // For many homogeneous ODEs, solution is implicit
        let f = expr!((v ^ 2) + 1);
        let solver = HomogeneousODESolver;
        let solution = solver.solve(&f, &y, &x);

        assert!(solution.is_ok());
        // Solution may be in implicit form F(x, v) = C
        let sol = solution.unwrap();
        let sol_str = sol.to_string();
        assert!(sol_str.contains("v") || sol_str.contains("x"));
    }
}