use crate::api::expr::Ex;
use crate::base::arena::Arena;
use crate::base::combinatorics::binomial;
use crate::base::node::{ExprId, ExprNode, SymbolId};
use crate::base::numeric::Q;
use crate::domains::matrix::Matrix;
use num_traits::One;
use num_traits::Signed;
const MAX_BRANCHES: usize = 8;
pub struct OdeResult {
pub solution: ExprId,
pub constants: Vec<ExprId>,
pub branches: Vec<ExprId>,
}
pub fn dsolve(
arena: &mut Arena,
expr: ExprId,
func: ExprId, var: ExprId, ) -> Option<OdeResult> {
let ExprNode::Symbol(func_sym) = *arena.node(func) else {
return None;
};
let first = dsolve_as_given(arena, expr, func, var);
if first
.as_ref()
.is_some_and(|r| !contains_sym(arena, r.solution, func_sym))
{
return first;
}
let mut best = first;
for form in normalized_forms(arena, expr, func, var, func_sym) {
let Some(result) = dsolve_as_given(arena, form, func, var) else {
continue;
};
if !contains_sym(arena, result.solution, func_sym) {
return Some(result);
}
if best.is_none() {
best = Some(result);
}
}
if let Some(r) = &best
&& let Some(&c) = r.constants.last()
&& let Some((solution, branches)) = explicit_families(arena, expr, r.solution, func, var, c)
{
return Some(OdeResult {
solution,
constants: r.constants.clone(),
branches,
});
}
best
}
fn normalized_forms(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Vec<ExprId> {
let mut forms: Vec<ExprId> = Vec::new();
let mut current = expr;
for _ in 0..4 {
match strip_nonderivative_factors(arena, current, func_sym) {
Some(next) if next != current => current = next,
_ => break,
}
}
if current != expr {
forms.push(current);
}
if let Some(lead) = divide_leading_coefficient(arena, current, func, var, func_sym)
&& lead != expr
&& !forms.contains(&lead)
{
forms.push(lead);
}
forms
}
fn has_derivative(arena: &Arena, e: ExprId) -> bool {
let mut stack = vec![e];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
let node = arena.node(n);
if matches!(node, ExprNode::Derivative(..)) {
return true;
}
node.for_each_child(|c| stack.push(c));
}
false
}
fn strip_nonderivative_factors(arena: &mut Arena, e: ExprId, func_sym: SymbolId) -> Option<ExprId> {
let positive_power_base = |arena: &Arena, f: ExprId| -> ExprId {
match arena.node(f) {
ExprNode::Pow(b, p) if arena.as_num(*p).is_some_and(|r| r.is_positive()) => *b,
_ => f,
}
};
match arena.node(e).clone() {
ExprNode::Mul(factors) => {
let carriers: Vec<ExprId> = factors
.iter()
.copied()
.filter(|&f| contains_sym(arena, f, func_sym) && has_derivative(arena, f))
.collect();
let [carrier] = carriers.as_slice() else {
return None;
};
let base = positive_power_base(arena, *carrier);
if base == *carrier && matches!(arena.node(*carrier), ExprNode::Pow(..)) {
return None;
}
Some(base)
}
ExprNode::Neg(inner) => Some(inner),
ExprNode::Pow(..) => {
let base = positive_power_base(arena, e);
(base != e && has_derivative(arena, base)).then_some(base)
}
_ => None,
}
}
fn divide_leading_coefficient(
arena: &mut Arena,
e: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<ExprId> {
const MAX_TERMS: usize = 64;
let chain = derivative_chain(arena, func, var);
let order = (1..chain.len())
.rev()
.find(|&k| expr_contains(arena, e, chain[k]))?;
let lead = chain[order];
let expanded = crate::transforms::expand::expand(arena, e);
let expanded = crate::transforms::eval::eval(arena, expanded);
let terms: Vec<ExprId> = match arena.node(expanded) {
ExprNode::Add(c) => c.to_vec(),
_ => return None,
};
if terms.len() > MAX_TERMS {
return None;
}
let mut groups: Vec<(ExprId, Vec<ExprId>)> = Vec::new();
for t in terms {
let (sign, inner) = match arena.node(t) {
ExprNode::Neg(i) => (true, *i),
_ => (false, t),
};
let factors: Vec<ExprId> = match arena.node(inner) {
ExprNode::Mul(fs) => fs.to_vec(),
_ => vec![inner],
};
let (in_func, coeff): (Vec<ExprId>, Vec<ExprId>) = factors
.into_iter()
.partition(|&f| contains_sym(arena, f, func_sym));
let key = match in_func.len() {
0 => arena.one,
1 => in_func[0],
_ => arena.mul(&in_func),
};
let mut c = match coeff.len() {
0 => arena.one,
1 => coeff[0],
_ => arena.mul(&coeff),
};
if sign {
c = arena.neg(c);
}
match groups.iter_mut().find(|(k, _)| *k == key) {
Some(g) => g.1.push(c),
None => groups.push((key, vec![c])),
}
}
let carriers: Vec<usize> = (0..groups.len())
.filter(|&i| expr_contains(arena, groups[i].0, lead))
.collect();
let [lead_group] = carriers.as_slice() else {
return None;
};
if groups[*lead_group].0 != lead {
return None;
}
let c = {
let parts = groups[*lead_group].1.clone();
let s = arena.add(&parts);
crate::transforms::eval::eval(arena, s)
};
if arena.as_num(c).is_some() || contains_sym(arena, c, func_sym) {
return None;
}
let inv = {
let m1 = arena.int(-1);
arena.pow(c, m1)
};
let mut out = Vec::with_capacity(groups.len());
for (key, parts) in groups {
let s = arena.add(&parts);
let q = arena.mul(&[s, inv]);
let q = crate::transforms::eval::eval(arena, q);
let q = crate::simplify::ratsimp::ratsimp(arena, q);
let q = crate::transforms::eval::eval(arena, q);
if arena.is_zero_structural(q) {
continue;
}
out.push(arena.mul(&[q, key]));
}
let sum = arena.add(&out);
Some(crate::transforms::eval::eval(arena, sum))
}
fn dsolve_as_given(
arena: &mut Arena,
expr: ExprId,
func: ExprId, var: ExprId, ) -> Option<OdeResult> {
let var_sym = match arena.node(var) {
ExprNode::Symbol(sid) => *sid,
_ => return None,
};
let func_sym = match arena.node(func) {
ExprNode::Symbol(sid) => *sid,
_ => return None,
};
let mut implicit: Option<OdeResult> = None;
macro_rules! attempt {
($result:expr) => {
if let Some(result) = $result {
if contains_sym(arena, result.solution, func_sym) {
if implicit.is_none() {
implicit = Some(result);
}
} else {
return Some(result);
}
}
};
}
attempt!(try_second_order_const_coeff(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_second_order_cc_nonhomogeneous(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_euler_cauchy(arena, expr, func, var, func_sym, var_sym));
attempt!(try_euler_cauchy_n(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_nth_order_linear_const_coeff(
arena, expr, func, var, func_sym
));
attempt!(try_linear_cc_symbolic(arena, expr, func, var, func_sym));
attempt!(try_variation_of_parameters(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_clairaut(arena, expr, func, var, func_sym, var_sym));
attempt!(try_first_order_linear_general(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_exact_ode(arena, expr, func, var, func_sym, var_sym));
attempt!(try_integrating_factor_ode(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_bernoulli(arena, expr, func, var, func_sym, var_sym));
attempt!(try_homogeneous_coefficient(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_nth_order_reducible(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_full_separable(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_simple_separable(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_rational_riccati(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_linear_argument(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_linear_coefficients(
arena, expr, func, var, func_sym, var_sym
));
attempt!(try_reduce_lowest_order(arena, expr, func, var, func_sym));
attempt!(try_solve_for_highest_derivative(
arena, expr, func, var, func_sym
));
if implicit.is_none() {
attempt!(try_first_order_linear_formal(
arena, expr, func, var, func_sym, var_sym
));
}
implicit
}
fn ode_order(arena: &mut Arena, expr: ExprId, func: ExprId, var: ExprId) -> Option<usize> {
let chain = derivative_chain(arena, func, var);
(1..chain.len())
.rev()
.find(|&k| expr_contains(arena, expr, chain[k]))
}
fn try_reduce_lowest_order(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<OdeResult> {
let n = ode_order(arena, expr, func, var)?;
if n < 3 {
return None;
}
let chain = derivative_chain(arena, func, var);
let mut probe = expr;
let mut placeholders = Vec::with_capacity(n);
for j in (1..=n).rev() {
let ph = arena.symbol(&format!("__red_d{j}"));
probe = crate::transforms::subs::subs(arena, probe, chain[j], ph);
placeholders.push((j, ph));
}
if contains_sym(arena, probe, func_sym) {
return None;
}
let k = placeholders
.iter()
.rev()
.find(|&&(_, ph)| crate::base::walk::contains(arena, probe, ph))
.map(|&(j, _)| j)?;
if k >= n {
return None;
}
let p = arena.symbol("__p_low");
let ExprNode::Symbol(p_sym) = *arena.node(p) else {
return None;
};
let pchain = derivative_chain(arena, p, var);
let mut reduced = expr;
for j in (k..=n).rev() {
reduced = crate::transforms::subs::subs(arena, reduced, chain[j], pchain[j - k]);
}
let reduced = crate::transforms::eval::eval(arena, reduced);
if contains_sym(arena, reduced, func_sym) {
return None;
}
let pr = dsolve(arena, reduced, p, var)?;
if contains_sym(arena, pr.solution, p_sym)
|| crate::base::walk::has_unevaluated(arena, pr.solution)
|| pr.constants.len() != n - k
{
return None;
}
let mut next = 1usize;
for &c in &pr.constants {
if let ExprNode::Symbol(s) = arena.node(c)
&& let Some(i) = arena
.symbol_name(*s)
.strip_prefix('C')
.and_then(|d| d.parse::<usize>().ok())
{
next = next.max(i + 1);
}
}
let fresh: Vec<ExprId> = (0..k)
.map(|j| arena.symbol(&format!("C{}", next + j)))
.collect();
let integrate_k = |arena: &mut Arena, fam: ExprId| -> Option<ExprId> {
let mut f = fam;
for step in 0..k {
let mut g = integrate_forms(arena, f, var);
if integral_failed(arena, g)
&& step > 0
&& let Some(h) = ln_abs_to_ln(arena, f)
{
g = integrate_forms(arena, h, var);
}
if integral_failed(arena, g) {
return None;
}
f = g;
}
let mut terms = vec![f];
for (j, &c) in fresh.iter().enumerate() {
let xj = match j {
0 => arena.one,
1 => var,
_ => {
let e = arena.int(i64::try_from(j).ok()?);
arena.pow(var, e)
}
};
terms.push(arena.mul(&[c, xj]));
}
let s = arena.add(&terms);
Some(crate::transforms::eval::eval(arena, s))
};
let solution = integrate_k(arena, pr.solution)?;
let branches = pr
.branches
.iter()
.filter_map(|&b| integrate_k(arena, b))
.collect();
let mut constants = pr.constants;
constants.extend(fresh);
Some(OdeResult {
solution,
constants,
branches,
})
}
const MAX_DERIVATIVE_ROOTS: usize = 4;
fn try_solve_for_highest_derivative(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<OdeResult> {
let n = ode_order(arena, expr, func, var)?;
let chain = derivative_chain(arena, func, var);
let t = arena.symbol("__t_hd");
let ExprNode::Symbol(t_sym) = *arena.node(t) else {
return None;
};
let g = crate::transforms::subs::subs(arena, expr, chain[n], t);
let g = crate::transforms::eval::eval(arena, g);
let coeffs = crate::transforms::solve::symbolic_poly_coeffs(arena, g, t)?;
if !(3..=MAX_DERIVATIVE_ROOTS + 1).contains(&coeffs.len()) {
return None;
}
let roots = crate::transforms::solve::solve(arena, g, t);
if roots.is_empty() || roots.iter().any(|r| contains_sym(arena, r.value, t_sym)) {
return None;
}
let mut families: Vec<ExprId> = Vec::new();
let mut constants: Option<Vec<ExprId>> = None;
for r in roots {
let neg = arena.neg(r.value);
let branch = arena.add(&[chain[n], neg]);
let branch = crate::transforms::eval::eval(arena, branch);
let Some(res) = dsolve(arena, branch, func, var) else {
continue;
};
if contains_sym(arena, res.solution, func_sym)
|| crate::base::walk::has_unevaluated(arena, res.solution)
|| res.constants.len() != n
{
continue;
}
match &constants {
Some(cs) if *cs != res.constants => continue,
Some(_) => {}
None => constants = Some(res.constants.clone()),
}
for f in std::iter::once(res.solution).chain(res.branches) {
if !families.contains(&f) && families.len() < MAX_BRANCHES {
families.push(f);
}
}
}
let constants = constants?;
if families.is_empty() {
return None;
}
let solution = families.remove(0);
Some(OdeResult {
solution,
constants,
branches: families,
})
}
fn first_order_rhs(arena: &mut Arena, expr: ExprId, func: ExprId, var: ExprId) -> Option<ExprId> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if expr_contains(arena, expr, d2y_dx2) {
return None;
}
let (m, n) = extract_m_n(arena, expr, func, var)?;
if expr_contains(arena, m, dy_dx) || expr_contains(arena, n, dy_dx) {
return None;
}
let n = crate::transforms::eval::eval(arena, n);
if arena.is_zero_structural(n) {
return None;
}
let neg_m = arena.neg(m);
let rhs = arena.div(neg_m, n);
Some(crate::transforms::eval::eval(arena, rhs))
}
fn affine_forms(
arena: &mut Arena,
e: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Vec<(Q, Q, ExprId)> {
use num_traits::Zero;
let mut out: Vec<(Q, Q, ExprId)> = Vec::new();
let post = crate::base::walk::post_order_ids(arena, e);
for id in post {
let ExprNode::Add(terms) = arena.node(id).clone() else {
continue;
};
let (mut a, mut b) = (Q::zero(), Q::zero());
let mut rest = Vec::new();
let mut affine = true;
for &t in terms.iter() {
let (c, term) = arena.as_coeff_term(t);
if term == var {
a += c;
} else if term == func {
b += c;
} else if contains_sym(arena, t, func_sym) || contains_sym(arena, t, var_sym) {
affine = false;
break;
} else {
rest.push(t);
}
}
if !affine || (a.is_zero() && b.is_zero()) {
continue;
}
let c = match rest.len() {
0 => arena.zero,
1 => rest[0],
_ => arena.add(&rest),
};
let c = crate::transforms::eval::eval(arena, c);
if !out
.iter()
.any(|(a2, b2, c2)| *a2 == a && *b2 == b && *c2 == c)
{
out.push((a, b, c));
}
}
out
}
fn try_linear_argument(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
use num_traits::Zero;
let rhs = first_order_rhs(arena, expr, func, var)?;
if !contains_sym(arena, rhs, func_sym) || !contains_sym(arena, rhs, var_sym) {
return None;
}
let u = arena.symbol("__u_lin");
let ExprNode::Symbol(u_sym) = *arena.node(u) else {
return None;
};
for (a, b, c) in affine_forms(arena, rhs, func, var, func_sym, var_sym) {
if a.is_zero() || b.is_zero() {
continue;
}
let a_e = ode_ratio_to_expr(arena, &a);
let b_e = ode_ratio_to_expr(arena, &b);
let ax = arena.mul(&[a_e, var]);
let neg_ax = arena.neg(ax);
let neg_c = arena.neg(c);
let num = arena.add(&[u, neg_ax, neg_c]);
let y_of_u = arena.div(num, b_e);
let y_of_u = crate::transforms::eval::eval(arena, y_of_u);
let g = crate::transforms::subs::subs(arena, rhs, func, y_of_u);
let mut g = crate::transforms::eval::eval(arena, g);
if contains_sym(arena, g, var_sym) {
g = crate::simplify::ratsimp::ratsimp(arena, g);
g = crate::transforms::eval::eval(arena, g);
}
if contains_sym(arena, g, var_sym) || !contains_sym(arena, g, u_sym) {
continue;
}
let du = arena.intern(ExprNode::Derivative(u, var));
let bg = arena.mul(&[b_e, g]);
let rate = arena.add(&[a_e, bg]);
let neg_rate = arena.neg(rate);
let u_ode = arena.add(&[du, neg_rate]);
let u_ode = crate::transforms::eval::eval(arena, u_ode);
let Some(ur) = dsolve(arena, u_ode, u, var) else {
continue;
};
if ur.constants.len() != 1 || crate::base::walk::has_unevaluated(arena, ur.solution) {
continue;
}
let by = arena.mul(&[b_e, func]);
let u_of_y = arena.add(&[ax, by, c]);
let back = |arena: &mut Arena, fam: ExprId| -> ExprId {
let r = if contains_sym(arena, fam, u_sym) {
crate::transforms::subs::subs(arena, fam, u, u_of_y)
} else {
crate::transforms::subs::subs(arena, y_of_u, u, fam)
};
crate::transforms::eval::eval(arena, r)
};
let solution = back(arena, ur.solution);
let branches = ur.branches.iter().map(|&f| back(arena, f)).collect();
return Some(OdeResult {
solution,
constants: ur.constants,
branches,
});
}
None
}
fn try_linear_coefficients(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
use num_traits::Zero;
let rhs = first_order_rhs(arena, expr, func, var)?;
if !contains_sym(arena, rhs, func_sym) || !contains_sym(arena, rhs, var_sym) {
return None;
}
let forms = affine_forms(arena, rhs, func, var, func_sym, var_sym);
let mut rational: Vec<(Q, Q, Q)> = Vec::with_capacity(forms.len());
for (a, b, c) in &forms {
rational.push((a.clone(), b.clone(), arena.as_num(*c)?.clone()));
}
if rational.iter().all(|(_, _, c)| c.is_zero()) {
return None;
}
let (i, j) = (0..rational.len())
.flat_map(|i| (i + 1..rational.len()).map(move |j| (i, j)))
.find(|&(i, j)| {
let (a1, b1, _) = &rational[i];
let (a2, b2, _) = &rational[j];
!(a1 * b2 - a2 * b1).is_zero()
})?;
let (a1, b1, c1) = &rational[i];
let (a2, b2, c2) = &rational[j];
let det = a1 * b2 - a2 * b1;
let h = (b1 * c2 - b2 * c1) / &det;
let k = (a2 * c1 - a1 * c2) / &det;
if h.is_zero() && k.is_zero() {
return None;
}
if rational
.iter()
.any(|(a, b, c)| !(a * &h + b * &k + c).is_zero())
{
return None;
}
let h_e = ode_ratio_to_expr(arena, &h);
let k_e = ode_ratio_to_expr(arena, &k);
let big_y = arena.symbol("__Y_lin");
let ExprNode::Symbol(big_y_sym) = *arena.node(big_y) else {
return None;
};
let x_shift = arena.add(&[var, h_e]);
let y_shift = arena.add(&[big_y, k_e]);
let shifted = crate::transforms::subs::subs(arena, rhs, func, y_shift);
let shifted = crate::transforms::subs::subs(arena, shifted, var, x_shift);
let shifted = crate::transforms::eval::eval(arena, shifted);
let shifted = crate::transforms::expand::expand(arena, shifted);
let shifted = crate::transforms::eval::eval(arena, shifted);
let dy = arena.intern(ExprNode::Derivative(big_y, var));
let neg = arena.neg(shifted);
let y_ode = arena.add(&[dy, neg]);
let y_ode = crate::transforms::eval::eval(arena, y_ode);
let yr = dsolve(arena, y_ode, big_y, var)?;
if yr.constants.len() != 1 || crate::base::walk::has_unevaluated(arena, yr.solution) {
return None;
}
let neg_h = arena.neg(h_e);
let x_back = arena.add(&[var, neg_h]);
let neg_k = arena.neg(k_e);
let y_back = arena.add(&[func, neg_k]);
let back = |arena: &mut Arena, fam: ExprId| -> ExprId {
let r = if contains_sym(arena, fam, big_y_sym) {
let r = crate::transforms::subs::subs(arena, fam, big_y, y_back);
crate::transforms::subs::subs(arena, r, var, x_back)
} else {
let r = crate::transforms::subs::subs(arena, fam, var, x_back);
arena.add(&[r, k_e])
};
crate::transforms::eval::eval(arena, r)
};
let solution = back(arena, yr.solution);
let branches = yr.branches.iter().map(|&f| back(arena, f)).collect();
Some(OdeResult {
solution,
constants: yr.constants,
branches,
})
}
fn try_simple_separable(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<OdeResult> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
if let ExprNode::Add(ref children) = arena.node(expr).clone() {
let mut deriv_coeffs: Vec<ExprId> = Vec::new();
let mut other_terms: Vec<ExprId> = Vec::new();
for &child in children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
deriv_coeffs.push(ode_ratio_to_expr(arena, &coeff));
} else if !contains_sym(arena, child, func_sym) {
other_terms.push(child);
} else if let ExprNode::Mul(factors) = arena.node(term).clone()
&& factors.iter().filter(|&&f| f == dy_dx).count() == 1
&& factors
.iter()
.all(|&f| f == dy_dx || !contains_sym(arena, f, func_sym))
{
let mut rest: Vec<ExprId> =
factors.iter().copied().filter(|&f| f != dy_dx).collect();
rest.push(ode_ratio_to_expr(arena, &coeff));
deriv_coeffs.push(arena.mul(&rest));
} else {
return None;
}
}
if !deriv_coeffs.is_empty() {
let a = arena.add(&deriv_coeffs);
let a = crate::transforms::eval::eval(arena, a);
if arena.is_zero_structural(a) {
return None;
}
let rhs = if other_terms.is_empty() {
arena.zero
} else {
let sum = arena.add(&other_terms);
let neg = arena.neg(sum);
let q = arena.div(neg, a);
crate::transforms::eval::eval(arena, q)
};
let integral = integrate_forms(arena, rhs, var);
let c1 = arena.symbol("C1");
let solution = arena.add(&[integral, c1]);
return Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
});
}
}
if expr == dy_dx {
let c1 = arena.symbol("C1");
return Some(OdeResult {
solution: c1,
constants: vec![c1],
branches: Vec::new(),
});
}
None
}
fn try_second_order_const_coeff(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
_func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<OdeResult> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut a_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut b_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut c_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == d2y_dx2 {
a_coeff += coeff;
} else if term == dy_dx {
b_coeff += coeff;
} else if term == func {
c_coeff += coeff;
} else {
return None; }
}
use num_traits::Zero;
if a_coeff.is_zero() {
return None; }
let b = &b_coeff / &a_coeff;
let c = &c_coeff / &a_coeff;
{
let four_r = num_rational::Ratio::<num_bigint::BigInt>::from_integer(4.into());
let disc = &b * &b - &four_r * &c;
if disc.is_negative() {
return build_trig_homogeneous_solution(arena, &b, &disc, var);
}
}
let r_var = arena.symbol("__r");
let two = arena.int(2);
let b_id = {
let nid = arena.intern_num(b.clone());
arena.intern(ExprNode::Num(nid))
};
let c_id = {
let nid = arena.intern_num(c.clone());
arena.intern(ExprNode::Num(nid))
};
let r_var_sq = arena.pow(r_var, two);
let b_r = arena.mul(&[b_id, r_var]);
let char_eq = arena.add(&[r_var_sq, b_r, c_id]);
let roots = crate::transforms::solve::solve(arena, char_eq, r_var);
let c1 = arena.symbol("C1");
let c2 = arena.symbol("C2");
match roots.len() {
2 => {
let r1 = roots[0].value;
let r2 = roots[1].value;
if r1 == r2 {
let rx = arena.mul(&[r1, var]);
let exp_rx = arena.exp(rx);
let c2_x = arena.mul(&[c2, var]);
let inner = arena.add(&[c1, c2_x]);
let solution = arena.mul(&[inner, exp_rx]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
} else {
let r1x = arena.mul(&[r1, var]);
let r2x = arena.mul(&[r2, var]);
let exp_r1x = arena.exp(r1x);
let exp_r2x = arena.exp(r2x);
let term1 = arena.mul(&[c1, exp_r1x]);
let term2 = arena.mul(&[c2, exp_r2x]);
let solution = arena.add(&[term1, term2]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
}
1 => {
let r = roots[0].value;
let rx = arena.mul(&[r, var]);
let exp_rx = arena.exp(rx);
let c2_x = arena.mul(&[c2, var]);
let inner = arena.add(&[c1, c2_x]);
let solution = arena.mul(&[inner, exp_rx]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
_ => None,
}
}
fn solve_characteristic_equation(arena: &mut Arena, b: Q, c: Q, var: ExprId) -> Option<OdeResult> {
{
let four_r = num_rational::Ratio::<num_bigint::BigInt>::from_integer(4.into());
let disc = &b * &b - &four_r * &c;
if disc.is_negative() {
return build_trig_homogeneous_solution(arena, &b, &disc, var);
}
}
let r_var = arena.symbol("__r");
let two = arena.int(2);
let b_id = {
let nid = arena.intern_num(b);
arena.intern(ExprNode::Num(nid))
};
let c_id = {
let nid = arena.intern_num(c);
arena.intern(ExprNode::Num(nid))
};
let r_var_sq = arena.pow(r_var, two);
let b_r = arena.mul(&[b_id, r_var]);
let char_eq = arena.add(&[r_var_sq, b_r, c_id]);
let roots = crate::transforms::solve::solve(arena, char_eq, r_var);
let c1 = arena.symbol("C1");
let c2 = arena.symbol("C2");
match roots.len() {
2 => {
let r1 = roots[0].value;
let r2 = roots[1].value;
if r1 == r2 {
let rx = arena.mul(&[r1, var]);
let exp_rx = arena.exp(rx);
let c2_x = arena.mul(&[c2, var]);
let inner = arena.add(&[c1, c2_x]);
let solution = arena.mul(&[inner, exp_rx]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
} else {
let r1x = arena.mul(&[r1, var]);
let r2x = arena.mul(&[r2, var]);
let exp_r1x = arena.exp(r1x);
let exp_r2x = arena.exp(r2x);
let term1 = arena.mul(&[c1, exp_r1x]);
let term2 = arena.mul(&[c2, exp_r2x]);
let solution = arena.add(&[term1, term2]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
}
1 => {
let r = roots[0].value;
let rx = arena.mul(&[r, var]);
let exp_rx = arena.exp(rx);
let c2_x = arena.mul(&[c2, var]);
let inner = arena.add(&[c1, c2_x]);
let solution = arena.mul(&[inner, exp_rx]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
_ => None,
}
}
fn try_second_order_cc_nonhomogeneous(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<OdeResult> {
use num_traits::Zero;
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut a_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut b_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut c_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut f_of_x_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == d2y_dx2 {
a_coeff += coeff;
} else if term == dy_dx {
b_coeff += coeff;
} else if term == func {
c_coeff += coeff;
} else if !contains_sym(arena, child, func_sym) {
f_of_x_terms.push(child);
} else {
return None; }
}
if a_coeff.is_zero() {
return None; }
if f_of_x_terms.is_empty() {
return None; }
let b = &b_coeff / &a_coeff;
let c = &c_coeff / &a_coeff;
let f_sum = if f_of_x_terms.len() == 1 {
f_of_x_terms[0]
} else {
arena.add(&f_of_x_terms)
};
let y_p_poly = if let Some(f_coeffs_expr) = arena.coefficients_of(f_sum, var) {
if f_coeffs_expr.is_empty() {
None
} else {
let mut rhs_coeffs: Vec<Q> = Vec::new();
let mut all_numeric = true;
for &cid in &f_coeffs_expr {
if let Some(val) = arena.as_num(cid) {
rhs_coeffs.push(-val.clone() / &a_coeff);
} else {
all_numeric = false;
break;
}
}
if all_numeric {
find_particular_polynomial(&b, &c, &rhs_coeffs)
.map(|pc| build_polynomial_expr(arena, &pc, var))
} else {
None
}
}
} else {
None
};
let y_p = if let Some(yp) = y_p_poly {
yp
} else {
let neg_f_sum = arena.neg(f_sum);
let a_id = ode_ratio_to_expr(arena, &a_coeff);
let rhs_expr = arena.div(neg_f_sum, a_id);
let rhs_expr = crate::transforms::eval::eval(arena, rhs_expr);
try_undetermined_trig_exp(arena, rhs_expr, &b, &c, var)?
};
let homo_result = solve_characteristic_equation(arena, b, c, var)?;
let solution = arena.add(&[homo_result.solution, y_p]);
Some(OdeResult {
solution,
constants: homo_result.constants,
branches: Vec::new(),
})
}
fn find_particular_polynomial(b: &Q, c: &Q, rhs_coeffs: &[Q]) -> Option<Vec<Q>> {
use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::Zero;
let n = rhs_coeffs.len() - 1;
if !c.is_zero() {
let mut a = vec![Ratio::<BigInt>::zero(); n + 1];
for j in (0..=n).rev() {
let a_j2 = if j + 2 <= n {
a[j + 2].clone()
} else {
Ratio::zero()
};
let a_j1 = if j < n {
a[j + 1].clone()
} else {
Ratio::zero()
};
let factor2 = Ratio::from_integer(BigInt::from(((j + 2) * (j + 1)) as i64)) * a_j2;
let factor1 = Ratio::from_integer(BigInt::from((j + 1) as i64)) * b.clone() * a_j1;
a[j] = (rhs_coeffs[j].clone() - factor2 - factor1) / c.clone();
}
Some(a)
} else if !b.is_zero() {
let mut bb = vec![Ratio::<BigInt>::zero(); n + 1];
for j in (0..=n).rev() {
let deriv_term = if j < n {
Ratio::from_integer(BigInt::from(((j + 2) * (j + 1)) as i64)) * bb[j + 1].clone()
} else {
Ratio::zero()
};
let denom = Ratio::from_integer(BigInt::from((j + 1) as i64)) * b.clone();
if denom.is_zero() {
return None;
}
bb[j] = (rhs_coeffs[j].clone() - deriv_term) / denom;
}
let mut result = vec![Ratio::zero()];
result.extend(bb);
Some(result)
} else {
let mut result = vec![Ratio::<BigInt>::zero(); 2];
for (j, r_j) in rhs_coeffs.iter().enumerate() {
let denom = Ratio::from_integer(BigInt::from(((j + 1) * (j + 2)) as i64));
result.push(r_j.clone() / denom);
}
Some(result)
}
}
fn build_polynomial_expr(arena: &mut Arena, coeffs: &[Q], var: ExprId) -> ExprId {
use num_traits::Zero;
let mut terms = Vec::new();
for (j, coeff) in coeffs.iter().enumerate() {
if coeff.is_zero() {
continue;
}
let x_pow_j = if j == 0 {
arena.one
} else if j == 1 {
var
} else {
let exp = arena.int(j as i64);
arena.pow(var, exp)
};
let term = arena.make_coeff_term(coeff.clone(), x_pow_j);
terms.push(term);
}
if terms.is_empty() {
arena.zero
} else if terms.len() == 1 {
terms[0]
} else {
arena.add(&terms)
}
}
fn try_first_order_linear(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<OdeResult> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut has_dy = false;
let mut a_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut f_of_x = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
has_dy = true;
if !coeff.is_one() {
return None; }
} else if term == func {
a_coeff += coeff;
} else if !contains_sym(arena, child, func_sym) {
f_of_x.push(child);
} else {
return None; }
}
use num_traits::Zero;
if !has_dy {
return None;
}
let a_id = {
let nid = arena.intern_num(a_coeff.clone());
arena.intern(ExprNode::Num(nid))
};
let neg_a = arena.neg(a_id);
let ax = arena.mul(&[a_id, var]);
let neg_ax = arena.mul(&[neg_a, var]);
let exp_ax = arena.exp(ax);
let exp_neg_ax = arena.exp(neg_ax);
if a_coeff.is_zero() {
let f = if f_of_x.is_empty() {
arena.zero
} else {
let sum = arena.add(&f_of_x);
arena.neg(sum)
};
let integral = crate::transforms::integrate::integrate(arena, f, var);
let c1 = arena.symbol("C1");
let solution = arena.add(&[integral, c1]);
return Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
});
}
let neg_f = if f_of_x.is_empty() {
arena.zero
} else {
let sum = arena.add(&f_of_x);
arena.neg(sum)
};
let integrand = if let ExprNode::Exp(neg_f_inner) = arena.node(neg_f).clone() {
let combined_arg = arena.add(&[neg_f_inner, ax]);
let combined_arg = crate::transforms::eval::eval(arena, combined_arg);
arena.exp(combined_arg)
} else {
arena.mul(&[neg_f, exp_ax])
};
let integrand = crate::transforms::eval::eval(arena, integrand);
let integral = crate::transforms::integrate::integrate(arena, integrand, var);
let c1 = arena.symbol("C1");
let inner = arena.add(&[integral, c1]);
let solution = arena.mul(&[exp_neg_ax, inner]);
Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
})
}
fn try_full_separable(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying full separable");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut has_dy = false;
let mut dy_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut other_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
has_dy = true;
dy_coeff += coeff;
} else {
other_terms.push(child);
}
}
use num_traits::Zero;
if !has_dy || dy_coeff.is_zero() {
return None;
}
if !dy_coeff.is_one() {
let inv_coeff = num_rational::Ratio::<num_bigint::BigInt>::one() / &dy_coeff;
let inv_id = {
let nid = arena.intern_num(inv_coeff);
arena.intern(ExprNode::Num(nid))
};
let mut scaled = Vec::new();
for &t in &other_terms {
scaled.push(arena.mul(&[inv_id, t]));
}
other_terms = scaled;
}
let rhs = if other_terms.is_empty() {
return None; } else {
let sum = arena.add(&other_terms);
arena.neg(sum)
};
if !contains_sym(arena, rhs, func_sym) {
return None; }
let mut factors = Vec::new();
for f in collect_mul_factors(arena, rhs) {
match arena.node(f).clone() {
ExprNode::Exp(arg) if matches!(arena.node(arg), ExprNode::Add(_)) => {
let ExprNode::Add(terms) = arena.node(arg).clone() else {
continue;
};
for &t in terms.iter() {
factors.push(arena.exp(t));
}
}
_ => factors.push(f),
}
}
let mut x_factors: Vec<ExprId> = Vec::new();
let mut y_factors: Vec<ExprId> = Vec::new();
for &factor in &factors {
let has_x = contains_sym(arena, factor, var_sym);
let has_y = contains_sym(arena, factor, func_sym);
if has_x && has_y {
return None;
} else if has_y {
y_factors.push(factor);
} else {
x_factors.push(factor);
}
}
if y_factors.is_empty() {
return None; }
let f_x = if x_factors.is_empty() {
arena.one
} else if x_factors.len() == 1 {
x_factors[0]
} else {
arena.mul(&x_factors)
};
let g_y = if y_factors.len() == 1 {
y_factors[0]
} else {
arena.mul(&y_factors)
};
{
let (coeff, base) = arena.as_coeff_term(g_y);
if base == func {
let coeff_id = ode_ratio_to_expr(arena, &coeff);
let scaled_fx = arena.mul(&[coeff_id, f_x]);
let integral_fx = crate::transforms::integrate::integrate(arena, scaled_fx, var);
let c1 = arena.symbol("C1");
let growth = arena.exp(integral_fx);
let solution = arena.mul(&[c1, growth]);
return Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
});
}
}
let neg_one = arena.int(-1);
let inv_gy = arena.pow(g_y, neg_one);
let inv_gy = crate::transforms::eval::eval(arena, inv_gy);
let lhs_attempt = integrate_forms(arena, inv_gy, func);
let lhs_integral = if integral_failed(arena, lhs_attempt) {
arena.intern(ExprNode::Integral(inv_gy, func))
} else {
lhs_attempt
};
let rhs_integral = crate::transforms::integrate::integrate(arena, f_x, var);
let c1 = arena.symbol("C1");
let neg_rhs = arena.neg(rhs_integral);
let neg_c1 = arena.neg(c1);
let implicit = arena.add(&[lhs_integral, neg_rhs, neg_c1]);
let implicit = crate::transforms::eval::eval(arena, implicit);
if !integral_failed(arena, implicit)
&& let Some((solution, branches)) = explicit_families(arena, expr, implicit, func, var, c1)
{
return Some(OdeResult {
solution,
constants: vec![c1],
branches,
});
}
Some(OdeResult {
solution: implicit,
constants: vec![c1],
branches: Vec::new(),
})
}
fn explicit_families(
arena: &mut Arena,
expr: ExprId,
implicit: ExprId,
func: ExprId,
var: ExprId,
c1: ExprId,
) -> Option<(ExprId, Vec<ExprId>)> {
let ExprNode::Symbol(func_sym) = *arena.node(func) else {
return None;
};
let explicit_roots = |arena: &mut Arena, e: ExprId| -> Option<Vec<ExprId>> {
let sols = crate::transforms::solve::solve(arena, e, func);
if sols.is_empty()
|| sols.len() > MAX_BRANCHES
|| sols.iter().any(|s| contains_sym(arena, s.value, func_sym))
{
return None;
}
let mut fams = Vec::with_capacity(sols.len());
for s in &sols {
let f = crate::transforms::eval::eval(arena, s.value);
if crate::base::walk::contains(arena, f, c1)
&& residual_verdict(arena, expr, f, func, var, c1) != Some(false)
{
fams.push(f);
}
}
(!fams.is_empty()).then_some(fams)
};
let stripped_families = |arena: &mut Arena| -> Option<Vec<ExprId>> {
let stripped = strip_log_abs(arena, implicit)?;
let fams = explicit_roots(arena, stripped)?;
let mut kept: Vec<ExprId> = Vec::with_capacity(fams.len());
for f in fams {
let f = absorb_exp_constant(arena, f, c1);
if !kept.contains(&f) && family_verifies(arena, expr, f, func, var, c1) {
kept.push(f);
}
}
(!kept.is_empty()).then_some(kept)
};
let mut fams = match explicit_roots(arena, implicit) {
Some(direct) if direct.iter().any(|&f| has_abs(arena, f)) => {
match stripped_families(arena) {
Some(s) if !s.iter().any(|&f| has_abs(arena, f)) => s,
_ => direct,
}
}
Some(direct) => direct,
None => stripped_families(arena)?,
};
for f in &mut fams {
let renamed = absorb_exp_constant(arena, *f, c1);
if renamed != *f && family_verifies(arena, expr, renamed, func, var, c1) {
*f = renamed;
}
*f = exp_log_to_power(arena, *f);
if let Some(renamed) = rename_constant_function(arena, *f, var, c1)
&& family_verifies(arena, expr, renamed, func, var, c1)
{
*f = renamed;
}
}
let neg_c1 = arena.neg(c1);
let mut unique: Vec<ExprId> = Vec::with_capacity(fams.len());
for f in fams {
let mirrored = crate::transforms::subs::subs(arena, f, c1, neg_c1);
let mirrored = crate::transforms::eval::eval(arena, mirrored);
if !unique.contains(&f) && !unique.contains(&mirrored) {
unique.push(f);
}
}
let first = unique.remove(0);
Some((first, unique))
}
fn exp_log_to_power(arena: &mut Arena, e: ExprId) -> ExprId {
let post = crate::base::walk::post_order_ids(arena, e);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
let mut changed = false;
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match *arena.node(rebuilt) {
ExprNode::Exp(arg) => {
let terms: Vec<ExprId> = match arena.node(arg) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![arg],
};
let mut powers = Vec::new();
let mut rest = Vec::new();
for t in terms {
let (k, term) = arena.as_coeff_term(t);
match *arena.node(term) {
ExprNode::Ln(u) => {
let ke = ode_ratio_to_expr(arena, &k);
powers.push(arena.pow(u, ke));
}
_ => rest.push(t),
}
}
if powers.is_empty() {
rebuilt
} else {
changed = true;
if !rest.is_empty() {
let r = arena.add(&rest);
powers.push(arena.exp(r));
}
arena.mul(&powers)
}
}
_ => rebuilt,
};
cache.insert(id, new);
}
if !changed {
return e;
}
let out = cache.get(&e).copied().unwrap_or(e);
crate::transforms::eval::eval(arena, out)
}
fn rename_constant_function(
arena: &mut Arena,
e: ExprId,
var: ExprId,
c1: ExprId,
) -> Option<ExprId> {
let ExprNode::Symbol(c_sym) = *arena.node(c1) else {
return None;
};
let ExprNode::Symbol(v_sym) = *arena.node(var) else {
return None;
};
let mut found: Option<ExprId> = None;
let mut stack = vec![e];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) || !contains_sym(arena, n, c_sym) {
continue;
}
if !contains_sym(arena, n, v_sym) {
match found {
Some(g) if g != n => return None,
_ => found = Some(n),
}
continue;
}
arena.node(n).for_each_child(|c| stack.push(c));
}
let g = found?;
if g == c1 || arena.as_num(g).is_some() {
return None;
}
if matches!(
arena.node(g),
ExprNode::Add(_) | ExprNode::Mul(_) | ExprNode::Neg(_)
) {
return None;
}
let out = crate::transforms::subs::subs(arena, e, g, c1);
Some(crate::transforms::eval::eval(arena, out))
}
fn has_abs(arena: &Arena, e: ExprId) -> bool {
let mut stack = vec![e];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
let node = arena.node(n);
if matches!(node, ExprNode::Abs(_)) {
return true;
}
node.for_each_child(|c| stack.push(c));
}
false
}
fn strip_log_abs(arena: &mut Arena, e: ExprId) -> Option<ExprId> {
let out = ln_abs_to_ln(arena, e)?;
let terms: Vec<ExprId> = match arena.node(out) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![out],
};
let mut lcm = num_bigint::BigInt::one();
for t in terms {
let (c, term) = arena.as_coeff_term(t);
if matches!(arena.node(term), ExprNode::Ln(_)) {
lcm = num_integer::Integer::lcm(&lcm, c.denom());
}
}
if lcm.is_one() {
return Some(out);
}
let scale = arena.intern_num(Q::from_integer(lcm));
let scale = arena.intern(ExprNode::Num(scale));
let scaled = arena.mul(&[scale, out]);
let scaled = crate::transforms::expand::expand(arena, scaled);
Some(crate::transforms::eval::eval(arena, scaled))
}
fn ln_abs_to_ln(arena: &mut Arena, e: ExprId) -> Option<ExprId> {
let post = crate::base::walk::post_order_ids(arena, e);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
let mut changed = false;
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match arena.node(rebuilt).clone() {
ExprNode::Ln(inner) => match arena.node(inner).clone() {
ExprNode::Abs(u) => {
changed = true;
arena.ln(u)
}
_ => rebuilt,
},
_ => rebuilt,
};
cache.insert(id, new);
}
if !changed {
return None;
}
let out = cache.get(&e).copied().unwrap_or(e);
Some(crate::transforms::eval::eval(arena, out))
}
fn absorb_exp_constant(arena: &mut Arena, e: ExprId, c1: ExprId) -> ExprId {
let ExprNode::Symbol(c_sym) = *arena.node(c1) else {
return e;
};
let split = |arena: &mut Arena, arg: ExprId| -> Option<(Q, Vec<ExprId>)> {
let terms: Vec<ExprId> = match arena.node(arg) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![arg],
};
let mut k: Option<Q> = None;
let mut rest = Vec::with_capacity(terms.len());
for t in terms {
let (c, term) = arena.as_coeff_term(t);
if term == c1 && k.is_none() {
k = Some(c);
} else if contains_sym(arena, t, c_sym) {
return None;
} else {
rest.push(t);
}
}
Some((k?, rest))
};
let post = crate::base::walk::post_order_ids(arena, e);
let mut erased: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
let mut hits: rustc_hash::FxHashMap<ExprId, (Q, Vec<ExprId>)> =
rustc_hash::FxHashMap::default();
let one = arena.one;
for &id in &post {
let gone = crate::base::walk::rebuild_with_cache(arena, id, &erased);
let hit = match *arena.node(id) {
ExprNode::Exp(arg) => split(arena, arg),
_ => None,
};
let new_gone = match hit {
Some(h) => {
hits.insert(id, h);
one
}
None => gone,
};
erased.insert(id, new_gone);
}
let rest = erased.get(&e).copied().unwrap_or(e);
if hits.is_empty() || contains_sym(arena, rest, c_sym) {
return e;
}
let mut num_gcd = num_bigint::BigInt::from(0);
let mut den_lcm = num_bigint::BigInt::one();
let mut sign_neg = None;
for &id in &post {
if let Some((k, _)) = hits.get(&id) {
num_gcd = num_integer::Integer::gcd(&num_gcd, k.numer());
den_lcm = num_integer::Integer::lcm(&den_lcm, k.denom());
sign_neg.get_or_insert(k.is_negative());
}
}
if num_gcd == num_bigint::BigInt::from(0) {
return e;
}
let mut k0 = Q::new(num_gcd, den_lcm);
if sign_neg == Some(true) {
k0 = -k0;
}
let mut renamed: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &renamed);
let new = match hits.get(&id).cloned() {
Some((k, rest)) => {
let r = arena.add(&rest);
let ex = arena.exp(r);
let p = ode_ratio_to_expr(arena, &(k / &k0));
let cp = arena.pow(c1, p);
arena.mul(&[cp, ex])
}
None => rebuilt,
};
renamed.insert(id, new);
}
let out = renamed.get(&e).copied().unwrap_or(e);
crate::transforms::eval::eval(arena, out)
}
fn family_verifies(
arena: &mut Arena,
expr: ExprId,
sol: ExprId,
func: ExprId,
var: ExprId,
c1: ExprId,
) -> bool {
match residual_verdict(arena, expr, sol, func, var, c1) {
Some(v) => v,
None => checkodesol(arena, expr, sol, func, var),
}
}
fn residual_verdict(
arena: &mut Arena,
expr: ExprId,
sol: ExprId,
func: ExprId,
var: ExprId,
c1: ExprId,
) -> Option<bool> {
let chain = derivative_chain(arena, func, var);
let order = (1..chain.len())
.rev()
.find(|&k| expr_contains(arena, expr, chain[k]))
.unwrap_or(0);
let mut ders = vec![sol];
for k in 1..=order {
let d = crate::transforms::diff::diff(arena, ders[k - 1], var);
ders.push(d);
}
let mut resid = expr;
for k in (0..=order).rev() {
resid = crate::transforms::subs::subs(arena, resid, chain[k], ders[k]);
}
let others: Vec<ExprId> = crate::base::walk::free_symbols(arena, sol)
.into_iter()
.filter(|&s| s != c1 && s != var && is_constant_name(arena, s))
.collect();
let (mut evaluated, mut passed) = (0usize, 0usize);
for (xn, xd) in [(3i64, 7i64), (11, 9), (-5, 4)] {
for (cn, cd) in [(2i64, 3i64), (-5, 4), (7, 5), (-1, 3)] {
let xv = arena.rational(xn, xd);
let cv = arena.rational(cn, cd);
let ov = arena.rational(3, 5);
let at = |arena: &mut Arena, e: ExprId| {
let mut e = crate::transforms::subs::subs(arena, e, var, xv);
e = crate::transforms::subs::subs(arena, e, c1, cv);
for &o in &others {
e = crate::transforms::subs::subs(arena, e, o, ov);
}
let e = crate::transforms::eval::eval(arena, e);
let z = crate::transforms::evalf::evalf_complex(arena, e, 12).ok()?;
crate::transforms::evalf::abs_to_f64(&z).filter(|m| m.is_finite())
};
let Some(r) = at(arena, resid) else {
continue;
};
let mut scale = 1.0f64;
for &d in &ders {
if let Some(m) = at(arena, d) {
scale = scale.max(m);
}
}
evaluated += 1;
if r <= 1e-9 * scale {
passed += 1;
if passed >= 2 {
return Some(true);
}
}
}
}
if passed >= 2 {
Some(true)
} else if evaluated >= 2 && passed == 0 {
Some(false)
} else {
None
}
}
fn is_constant_name(arena: &Arena, s: ExprId) -> bool {
match arena.node(s) {
ExprNode::Symbol(sid) => {
let name = arena.symbol_name(*sid);
name.len() > 1 && name.starts_with('C') && name[1..].chars().all(|c| c.is_ascii_digit())
}
_ => false,
}
}
fn collect_mul_factors(arena: &mut Arena, expr: ExprId) -> Vec<ExprId> {
let mut factors = Vec::new();
let mut expr = expr;
loop {
match arena.node(expr).clone() {
ExprNode::Mul(children) => {
factors.extend(children.iter().copied());
return factors;
}
ExprNode::Neg(inner) => {
factors.push(arena.int(-1));
expr = inner;
}
_ => {
factors.push(expr);
return factors;
}
}
}
}
fn try_first_order_linear_general(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
first_order_linear(arena, expr, func, var, func_sym, var_sym, false)
}
fn try_first_order_linear_formal(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
first_order_linear(arena, expr, func, var, func_sym, var_sym, true)
}
fn first_order_linear(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
formal: bool,
) -> Option<OdeResult> {
tracing::debug!("ode: trying variable-coefficient first-order linear");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut has_dy = false;
let mut dy_coeff_rational = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut y_terms: Vec<ExprId> = Vec::new(); let mut free_terms: Vec<ExprId> = Vec::new(); let mut dy_coeff_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
has_dy = true;
dy_coeff_rational += coeff;
} else if expr_contains(arena, child, dy_dx) {
let ExprNode::Mul(factors) = arena.node(term).clone() else {
return None;
};
let (with_dy, rest): (Vec<ExprId>, Vec<ExprId>) =
factors.iter().copied().partition(|&f| f == dy_dx);
if with_dy.len() != 1 || rest.iter().any(|&f| contains_sym(arena, f, func_sym)) {
return None;
}
has_dy = true;
let c = ode_ratio_to_expr(arena, &coeff);
let mut parts = vec![c];
parts.extend(rest);
dy_coeff_terms.push(arena.mul(&parts));
} else if contains_sym(arena, child, func_sym) {
y_terms.push(child);
} else {
free_terms.push(child);
}
}
use num_traits::Zero;
if !has_dy {
return None;
}
let dy_coeff_expr: Option<ExprId> = if dy_coeff_terms.is_empty() {
None
} else {
let mut parts = dy_coeff_terms.clone();
if !dy_coeff_rational.is_zero() {
parts.push(ode_ratio_to_expr(arena, &dy_coeff_rational));
}
let a = arena.add(&parts);
let a = crate::transforms::eval::eval(arena, a);
match arena.as_num(a) {
Some(r) => {
dy_coeff_rational = r.clone();
None
}
None => Some(a),
}
};
if dy_coeff_expr.is_none() && dy_coeff_rational.is_zero() {
return None;
}
if y_terms.is_empty() {
return None;
}
let mut p_x_terms: Vec<ExprId> = Vec::new();
for &yt in &y_terms {
{
let px = extract_coeff_of_func(arena, yt, func, func_sym, var_sym)?;
p_x_terms.push(px);
}
}
if let Some(a) = dy_coeff_expr {
let neg_one = arena.int(-1);
let inv_a = arena.pow(a, neg_one);
p_x_terms = p_x_terms
.iter()
.map(|&p| {
let t = arena.mul(&[inv_a, p]);
crate::transforms::eval::eval(arena, t)
})
.collect();
free_terms = free_terms
.iter()
.map(|&f| {
let t = arena.mul(&[inv_a, f]);
crate::transforms::eval::eval(arena, t)
})
.collect();
} else if !dy_coeff_rational.is_one() {
let inv_coeff = num_rational::Ratio::<num_bigint::BigInt>::one() / &dy_coeff_rational;
let inv_id = {
let nid = arena.intern_num(inv_coeff);
arena.intern(ExprNode::Num(nid))
};
let mut scaled_p = Vec::new();
for &p in &p_x_terms {
scaled_p.push(arena.mul(&[inv_id, p]));
}
p_x_terms = scaled_p;
let mut scaled_f = Vec::new();
for &f in &free_terms {
scaled_f.push(arena.mul(&[inv_id, f]));
}
free_terms = scaled_f;
}
let p_x = if p_x_terms.is_empty() {
arena.zero
} else if p_x_terms.len() == 1 {
p_x_terms[0]
} else {
arena.add(&p_x_terms)
};
if contains_sym(arena, p_x, func_sym) {
return None;
}
let p_x = crate::transforms::eval::eval(arena, p_x);
if arena.as_num(p_x).is_some() {
let normalised = if dy_coeff_expr.is_none() && dy_coeff_rational.is_one() {
expr
} else {
let py = arena.mul(&[p_x, func]);
let mut parts = vec![dy_dx, py];
parts.extend(free_terms.iter().copied());
let e = arena.add(&parts);
crate::transforms::eval::eval(arena, e)
};
if let Some(r) = try_first_order_linear(arena, normalised, func, var, func_sym, var_sym) {
return Some(r);
}
}
let q_x = if free_terms.is_empty() {
arena.zero
} else {
let sum = arena.add(&free_terms);
arena.neg(sum)
};
let int_px = crate::transforms::integrate::integrate(arena, p_x, var);
let mu = if integral_failed(arena, int_px) {
if !formal {
return None;
}
let p_x = crate::transforms::eval::eval(arena, p_x);
let formal_int = arena.intern(ExprNode::Integral(p_x, var));
arena.exp(formal_int)
} else if formal {
return None; } else {
exp_of_log_sum(arena, int_px)
};
let neg_one = arena.int(-1);
let inv_mu = arena.pow(mu, neg_one);
let c1 = arena.symbol("C1");
if arena.is_zero_structural(q_x) {
let solution = arena.mul(&[c1, inv_mu]);
return Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
});
}
let integrand = arena.mul(&[q_x, mu]);
let integrand = crate::transforms::eval::eval(arena, integrand);
let integral = integrate_forms(arena, integrand, var);
let inner = arena.add(&[integral, c1]);
let solution = arena.mul(&[inv_mu, inner]);
Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
})
}
fn extract_coeff_of_func(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<ExprId> {
if expr == func {
return Some(arena.one);
}
if let ExprNode::Neg(inner) = arena.node(expr).clone() {
if inner == func {
let neg_one = arena.int(-1);
return Some(neg_one);
}
if let Some(inner_coeff) = extract_coeff_of_func(arena, inner, func, func_sym, _var_sym) {
let result = arena.neg(inner_coeff);
return Some(result);
}
return None;
}
if let ExprNode::Mul(ref children) = arena.node(expr).clone() {
let mut found_y = false;
let mut other_factors: Vec<ExprId> = Vec::new();
let mut y_count = 0;
for &child in children {
if child == func {
y_count += 1;
if y_count > 1 {
return None; }
found_y = true;
} else if contains_sym(arena, child, func_sym) {
return None;
} else {
other_factors.push(child);
}
}
if found_y {
let coeff = if other_factors.is_empty() {
arena.one
} else if other_factors.len() == 1 {
other_factors[0]
} else {
arena.mul(&other_factors)
};
return Some(coeff);
}
}
{
let (coeff, term) = arena.as_coeff_term(expr);
if !coeff.is_one()
&& term != expr
&& let Some(inner_coeff) = extract_coeff_of_func(arena, term, func, func_sym, _var_sym)
{
let coeff_id = {
let nid = arena.intern_num(coeff);
arena.intern(ExprNode::Num(nid))
};
let result = arena.mul(&[coeff_id, inner_coeff]);
return Some(result);
}
}
None
}
fn extract_m_n(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
) -> Option<(ExprId, ExprId)> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
if expr == dy_dx {
return Some((arena.zero, arena.one));
}
let children = match arena.node(expr).clone() {
ExprNode::Add(c) => c,
_ => return None,
};
let mut m_terms: Vec<ExprId> = Vec::new();
let mut n_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
let coeff_id = {
let nid = arena.intern_num(coeff);
arena.intern(ExprNode::Num(nid))
};
n_terms.push(coeff_id);
} else if expr_contains(arena, child, dy_dx) {
if let ExprNode::Mul(ref mul_children) = arena.node(child).clone() {
let mut found_dy = false;
let mut other_factors: Vec<ExprId> = Vec::new();
for &mc in mul_children.iter() {
if mc == dy_dx && !found_dy {
found_dy = true;
} else {
other_factors.push(mc);
}
}
if found_dy {
let n_factor = match other_factors.len() {
0 => arena.one,
1 => other_factors[0],
_ => arena.mul(&other_factors),
};
n_terms.push(n_factor);
} else {
return None; }
} else {
return None;
}
} else {
m_terms.push(child);
}
}
if n_terms.is_empty() {
return None; }
let m_expr = match m_terms.len() {
0 => arena.zero,
1 => m_terms[0],
_ => arena.add(&m_terms),
};
let n_expr = match n_terms.len() {
1 => n_terms[0],
_ => arena.add(&n_terms),
};
Some((m_expr, n_expr))
}
fn try_exact_ode(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying exact ODE");
let (m_expr, n_expr) = extract_m_n(arena, expr, func, var)?;
if !contains_sym(arena, m_expr, func_sym) && !contains_sym(arena, n_expr, func_sym) {
return None;
}
let dm_dy = crate::transforms::diff::diff(arena, m_expr, func);
let dn_dx = crate::transforms::diff::diff(arena, n_expr, var);
let diff_check = arena.sub(dm_dy, dn_dx);
let diff_eval = crate::transforms::eval::eval(arena, diff_check);
let diff_expanded = crate::transforms::expand::expand(arena, diff_eval);
let diff_simplified = crate::transforms::eval::eval(arena, diff_expanded);
if diff_simplified != arena.zero && !rational_zero(arena, diff_simplified) {
return None; }
let integral_m = crate::transforms::integrate::integrate(arena, m_expr, var);
if integral_failed(arena, integral_m) {
return None; }
let integral_m = generic_branches(arena, integral_m, var).unwrap_or(integral_m);
let g_prime_of = |arena: &mut Arena, im: ExprId| -> ExprId {
let d_intm_dy = crate::transforms::diff::diff(arena, im, func);
let g_prime = arena.sub(n_expr, d_intm_dy);
let g_prime = crate::transforms::eval::eval(arena, g_prime);
let g_prime = crate::transforms::expand::expand(arena, g_prime);
crate::transforms::eval::eval(arena, g_prime)
};
let mut integral_m = integral_m;
let mut g_prime = g_prime_of(arena, integral_m);
if contains_sym(arena, g_prime, var_sym)
&& let Some(stripped) = ln_abs_to_ln(arena, integral_m)
{
let g2 = g_prime_of(arena, stripped);
if !contains_sym(arena, g2, var_sym) {
integral_m = stripped;
g_prime = g2;
}
}
if contains_sym(arena, g_prime, var_sym) {
return None;
}
let g_y = crate::transforms::integrate::integrate(arena, g_prime, func);
if integral_failed(arena, g_y) {
return None;
}
let potential = arena.add(&[integral_m, g_y]);
let potential = crate::transforms::eval::eval(arena, potential);
let c1 = arena.symbol("C1");
let f_minus_c1 = arena.sub(potential, c1);
let solutions = crate::transforms::solve::solve(arena, f_minus_c1, func);
if solutions.len() == 1 {
return Some(OdeResult {
solution: solutions[0].value,
constants: vec![c1],
branches: Vec::new(),
});
}
let implicit = crate::transforms::eval::eval(arena, f_minus_c1);
Some(OdeResult {
solution: implicit,
constants: vec![c1],
branches: Vec::new(),
})
}
fn try_integrating_factor_ode(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying integrating factor for non-exact ODE");
let (m_expr, n_expr) = extract_m_n(arena, expr, func, var)?;
if !contains_sym(arena, m_expr, func_sym) && !contains_sym(arena, n_expr, func_sym) {
return None;
}
let dm_dy = crate::transforms::diff::diff(arena, m_expr, func);
let dn_dx = crate::transforms::diff::diff(arena, n_expr, var);
let diff_mn = arena.sub(dm_dy, dn_dx); let diff_eval = crate::transforms::eval::eval(arena, diff_mn);
let diff_expanded = crate::transforms::expand::expand(arena, diff_eval);
let diff_simplified = crate::transforms::eval::eval(arena, diff_expanded);
if diff_simplified == arena.zero || rational_zero(arena, diff_simplified) {
return try_exact_ode(arena, expr, func, var, func_sym, var_sym);
}
{
let ratio = arena.div(diff_simplified, n_expr);
let ratio = crate::transforms::eval::eval(arena, ratio);
let ratio = crate::transforms::expand::expand(arena, ratio);
let ratio = crate::transforms::eval::eval(arena, ratio);
let ratio_cancelled = crate::simplify::ratsimp::ratsimp(arena, ratio);
let ratio_cancelled = crate::transforms::eval::eval(arena, ratio_cancelled);
if !contains_sym(arena, ratio_cancelled, func_sym) {
let int_ratio = crate::transforms::integrate::integrate(arena, ratio_cancelled, var);
if !integral_failed(arena, int_ratio) {
let mu = exp_of_log_sum(arena, int_ratio);
let new_m = arena.mul(&[mu, m_expr]);
let new_n = arena.mul(&[mu, n_expr]);
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let n_dy = arena.mul(&[new_n, dy_dx]);
let new_expr = arena.add(&[new_m, n_dy]);
let new_expr = crate::transforms::eval::eval(arena, new_expr);
if let Some(result) = try_exact_ode(arena, new_expr, func, var, func_sym, var_sym) {
return Some(result);
}
}
}
}
{
let neg_diff = arena.neg(diff_simplified); let ratio = arena.div(neg_diff, m_expr);
let ratio = crate::transforms::eval::eval(arena, ratio);
let ratio = crate::transforms::expand::expand(arena, ratio);
let ratio = crate::transforms::eval::eval(arena, ratio);
let ratio_cancelled = crate::simplify::ratsimp::ratsimp(arena, ratio);
let ratio_cancelled = crate::transforms::eval::eval(arena, ratio_cancelled);
if !contains_sym(arena, ratio_cancelled, var_sym) {
let int_ratio = crate::transforms::integrate::integrate(arena, ratio_cancelled, func);
if !integral_failed(arena, int_ratio) {
let mu = exp_of_log_sum(arena, int_ratio);
let new_m = arena.mul(&[mu, m_expr]);
let new_n = arena.mul(&[mu, n_expr]);
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let n_dy = arena.mul(&[new_n, dy_dx]);
let new_expr = arena.add(&[new_m, n_dy]);
let new_expr = crate::transforms::eval::eval(arena, new_expr);
if let Some(result) = try_exact_ode(arena, new_expr, func, var, func_sym, var_sym) {
return Some(result);
}
}
}
}
None
}
fn try_homogeneous_coefficient(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying homogeneous coefficient");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if expr_contains(arena, expr, d2y_dx2) {
return None;
}
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut has_dy = false;
let mut dy_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut other_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
has_dy = true;
dy_coeff += coeff;
} else {
other_terms.push(child);
}
}
use num_traits::Zero;
if !has_dy || dy_coeff.is_zero() {
return None;
}
let rhs = if other_terms.is_empty() {
return None;
} else {
let sum = arena.add(&other_terms);
arena.neg(sum)
};
let rhs = if !dy_coeff.is_one() {
let inv_id = ode_ratio_to_expr(
arena,
&(num_rational::Ratio::<num_bigint::BigInt>::one() / &dy_coeff),
);
let s = arena.mul(&[inv_id, rhs]);
crate::transforms::eval::eval(arena, s)
} else {
rhs
};
if !contains_sym(arena, rhs, func_sym) || !contains_sym(arena, rhs, var_sym) {
return None;
}
let v = arena.symbol("__v");
let rhs_sub = crate::transforms::subs::subs(arena, rhs, func, v);
let rhs_sub = crate::transforms::subs::subs(arena, rhs_sub, var, arena.one);
let rhs_sub = crate::transforms::eval::eval(arena, rhs_sub);
let rhs_sub = crate::transforms::expand::expand(arena, rhs_sub);
let rhs_sub = crate::transforms::eval::eval(arena, rhs_sub);
if contains_sym(arena, rhs_sub, var_sym) {
return None;
}
{
let vx = arena.mul(&[v, var]);
let probe = crate::transforms::subs::subs(arena, rhs, func, vx);
let probe = crate::transforms::eval::eval(arena, probe);
let diff = arena.sub(probe, rhs_sub);
let diff = crate::transforms::eval::eval(arena, diff);
let diff = crate::simplify::ratsimp::ratsimp(arena, diff);
let diff = crate::transforms::eval::eval(arena, diff);
if !arena.is_zero_structural(diff) {
return None;
}
}
let f_v_minus_v = arena.sub(rhs_sub, v);
let f_v_minus_v = crate::transforms::eval::eval(arena, f_v_minus_v);
let f_v_minus_v = crate::simplify::ratsimp::ratsimp(arena, f_v_minus_v);
let f_v_minus_v = crate::transforms::eval::eval(arena, f_v_minus_v);
if f_v_minus_v == arena.zero {
let c1 = arena.symbol("C1");
let solution = arena.mul(&[c1, var]);
return Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
});
}
let neg_one = arena.int(-1);
let inv_fv = arena.pow(f_v_minus_v, neg_one);
let lhs_integral = integrate_forms(arena, inv_fv, v);
if integral_failed(arena, lhs_integral) {
return None;
}
let abs_x = arena.abs(var);
let ln_abs_x = arena.ln(abs_x);
let c1 = arena.symbol("C1");
let rhs_eq = arena.add(&[ln_abs_x, c1]);
let y_over_x = arena.div(func, var);
let lhs_backsub = crate::transforms::subs::subs(arena, lhs_integral, v, y_over_x);
let lhs_backsub = crate::transforms::eval::eval(arena, lhs_backsub);
let implicit = arena.sub(lhs_backsub, rhs_eq);
let implicit = crate::transforms::eval::eval(arena, implicit);
let solutions = crate::transforms::solve::solve(arena, implicit, func);
if solutions.len() == 1 {
let sol = crate::transforms::eval::eval(arena, solutions[0].value);
return Some(OdeResult {
solution: sol,
constants: vec![c1],
branches: Vec::new(),
});
}
Some(OdeResult {
solution: implicit,
constants: vec![c1],
branches: Vec::new(),
})
}
fn try_nth_order_reducible(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
_func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying nth-order reducible (missing x)");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if !expr_contains(arena, expr, d2y_dx2) {
return None;
}
let d2_placeholder = arena.symbol("__d2");
let d1_placeholder = arena.symbol("__d1");
let stripped = crate::transforms::subs::subs(arena, expr, d2y_dx2, d2_placeholder);
let stripped = crate::transforms::subs::subs(arena, stripped, dy_dx, d1_placeholder);
let has_x = contains_sym(arena, stripped, var_sym);
let has_y = contains_sym(arena, stripped, _func_sym);
if has_x && has_y {
return None; }
let p = arena.symbol("__p");
let p_sym = match arena.node(p) {
ExprNode::Symbol(s) => *s,
_ => return None,
};
let c2 = arena.symbol("C2");
let usable = |arena: &Arena, r: &OdeResult| {
!contains_sym(arena, r.solution, p_sym)
&& r.constants.len() == 1
&& expr_contains(arena, r.solution, r.constants[0])
&& !crate::base::walk::has_unevaluated(arena, r.solution)
};
if !has_y {
let dp_dx = arena.intern(ExprNode::Derivative(p, var));
let reduced = crate::transforms::subs::subs(arena, expr, d2y_dx2, dp_dx);
let reduced = crate::transforms::subs::subs(arena, reduced, dy_dx, p);
let reduced = crate::transforms::eval::eval(arena, reduced);
if let Some(pr) = dsolve(arena, reduced, p, var)
&& usable(arena, &pr)
{
let antiderivative = |arena: &mut Arena, p_family: ExprId| {
let integral = crate::transforms::integrate::integrate(arena, p_family, var);
(!integral_failed(arena, integral)).then(|| {
let s = arena.add(&[integral, c2]);
crate::transforms::eval::eval(arena, s)
})
};
if let Some(solution) = antiderivative(arena, pr.solution) {
let mut branches = Vec::new();
for &b in &pr.branches {
match antiderivative(arena, b) {
Some(y_b) => branches.push(y_b),
None => tracing::debug!("ode: a family of y' has no closed antiderivative"),
}
}
return Some(OdeResult {
solution,
constants: vec![pr.constants[0], c2],
branches,
});
}
}
if has_x {
return None;
}
}
let dp_dy = arena.intern(ExprNode::Derivative(p, func));
let p_dp_dy = arena.mul(&[p, dp_dy]);
let reduced = crate::transforms::subs::subs(arena, expr, d2y_dx2, p_dp_dy);
let reduced = crate::transforms::subs::subs(arena, reduced, dy_dx, p);
let reduced = crate::transforms::eval::eval(arena, reduced);
let reduced = divide_common_factor(arena, reduced, p).unwrap_or(reduced);
let p_result = dsolve(arena, reduced, p, func)?;
if !usable(arena, &p_result) {
return None;
}
let p_sol = p_result.solution;
let neg_one_id = arena.int(-1);
let inv_p = arena.pow(p_sol, neg_one_id);
let inv_p = crate::transforms::eval::eval(arena, inv_p);
let lhs_integral = crate::transforms::integrate::integrate(arena, inv_p, func);
if integral_failed(arena, lhs_integral) {
return None; }
let rhs = arena.add(&[var, c2]);
let implicit = arena.sub(lhs_integral, rhs);
let implicit = crate::transforms::eval::eval(arena, implicit);
let solutions = crate::transforms::solve::solve(arena, implicit, func);
if solutions.len() == 1 {
let sol = crate::transforms::eval::eval(arena, solutions[0].value);
let mut constants = p_result.constants;
constants.push(c2);
return Some(OdeResult {
solution: sol,
constants,
branches: Vec::new(),
});
}
let mut constants = p_result.constants;
constants.push(c2);
Some(OdeResult {
solution: implicit,
constants,
branches: Vec::new(),
})
}
fn contains_sym(arena: &Arena, expr: ExprId, sym: SymbolId) -> bool {
crate::base::walk::has_free_symbol(arena, expr, sym)
}
fn integrate_forms(arena: &mut Arena, e: ExprId, var: ExprId) -> ExprId {
let first = crate::transforms::integrate::integrate(arena, e, var);
if !integral_failed(arena, first) {
return first;
}
let merged = crate::simplify::powsimp::powsimp(arena, e);
let merged = crate::transforms::eval::eval(arena, merged);
let fraction = crate::simplify::ratsimp::ratsimp(arena, merged);
let fraction = crate::transforms::eval::eval(arena, fraction);
let expanded = crate::transforms::expand::expand(arena, merged);
let expanded = crate::transforms::eval::eval(arena, expanded);
let factored = {
let post = crate::base::walk::post_order_ids(arena, merged);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = if matches!(arena.node(rebuilt), ExprNode::Add(_)) {
crate::simplify::factor_terms::factor_terms(arena, rebuilt)
} else {
rebuilt
};
cache.insert(id, new);
}
let f = cache.get(&merged).copied().unwrap_or(merged);
crate::transforms::eval::eval(arena, f)
};
for form in [merged, fraction, expanded, factored] {
if form == e {
continue;
}
let r = crate::transforms::integrate::integrate(arena, form, var);
if !integral_failed(arena, r) {
return r;
}
}
first
}
fn rational_zero(arena: &mut Arena, e: ExprId) -> bool {
let r = crate::simplify::ratsimp::ratsimp(arena, e);
let r = crate::transforms::eval::eval(arena, r);
arena.is_zero_structural(r)
}
fn divide_common_factor(arena: &mut Arena, expr: ExprId, p: ExprId) -> Option<ExprId> {
let expanded = crate::transforms::expand::expand(arena, expr);
let terms: Vec<ExprId> = match arena.node(expanded) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![expanded],
};
let has_factor = |arena: &Arena, t: ExprId| -> bool {
let is_power_of_p = |f: ExprId| -> bool {
f == p
|| matches!(arena.node(f), ExprNode::Pow(b, e) if *b == p
&& arena.as_num(*e).is_some_and(|r| r.is_integer() && r.is_positive()))
};
match arena.node(t) {
ExprNode::Mul(fs) => fs.iter().any(|&f| is_power_of_p(f)),
ExprNode::Neg(inner) => match arena.node(*inner) {
ExprNode::Mul(fs) => fs.iter().any(|&f| is_power_of_p(f)),
_ => is_power_of_p(*inner),
},
_ => is_power_of_p(t),
}
};
if terms.is_empty() || !terms.iter().all(|&t| has_factor(arena, t)) {
return None;
}
let quotient = arena.div(expanded, p);
let quotient = crate::transforms::eval::eval(arena, quotient);
let quotient = crate::transforms::expand::expand(arena, quotient);
Some(crate::transforms::eval::eval(arena, quotient))
}
fn integral_failed(arena: &Arena, id: ExprId) -> bool {
let mut stack = vec![id];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
let node = arena.node(n);
if matches!(node, ExprNode::Integral(..)) {
return true;
}
node.for_each_child(|c| stack.push(c));
}
false
}
fn exp_of_log_sum(arena: &mut Arena, integral: ExprId) -> ExprId {
let terms: Vec<ExprId> = match arena.node(integral).clone() {
ExprNode::Add(c) => c.to_vec(),
_ => vec![integral],
};
let mut factors: Vec<ExprId> = Vec::new();
let mut leftover: Vec<ExprId> = Vec::new();
for t in terms {
let (coeff, term) = arena.as_coeff_term(t);
let inner = match arena.node(term).clone() {
ExprNode::Ln(inner) => Some(inner),
_ => None,
};
match inner {
Some(inner) => {
let base = match arena.node(inner).clone() {
ExprNode::Abs(a) => a,
_ => inner,
};
if coeff.is_one() {
factors.push(base);
} else {
let c_id = ode_ratio_to_expr(arena, &coeff);
factors.push(arena.pow(base, c_id));
}
}
None => leftover.push(t),
}
}
if factors.is_empty() {
return arena.exp(integral);
}
if !leftover.is_empty() {
let rest = arena.add(&leftover);
factors.push(arena.exp(rest));
}
let prod = arena.mul(&factors);
crate::transforms::eval::eval(arena, prod)
}
fn ode_ratio_to_expr(arena: &mut Arena, r: &Q) -> ExprId {
let nid = arena.intern_num(r.clone());
arena.intern(ExprNode::Num(nid))
}
fn build_trig_homogeneous_solution(
arena: &mut Arena,
b: &Q,
disc: &Q,
var: ExprId,
) -> Option<OdeResult> {
use num_traits::Zero;
let c1 = arena.symbol("C1");
let c2 = arena.symbol("C2");
let two_r = num_rational::Ratio::<num_bigint::BigInt>::from_integer(2.into());
let alpha = -(b.clone()) / &two_r;
let neg_disc = -(disc.clone());
let neg_disc_id = ode_ratio_to_expr(arena, &neg_disc);
let half = arena.rational(1, 2);
let sqrt_neg_disc = arena.pow(neg_disc_id, half);
let two_id = arena.int(2);
let beta_expr = arena.div(sqrt_neg_disc, two_id);
let beta_expr = crate::transforms::eval::eval(arena, beta_expr);
let beta_x = arena.mul(&[beta_expr, var]);
let cos_bx = arena.cos(beta_x);
let sin_bx = arena.sin(beta_x);
let c1_cos = arena.mul(&[c1, cos_bx]);
let c2_sin = arena.mul(&[c2, sin_bx]);
let trig_part = arena.add(&[c1_cos, c2_sin]);
let solution = if alpha.is_zero() {
trig_part
} else {
let alpha_id = ode_ratio_to_expr(arena, &alpha);
let alpha_x = arena.mul(&[alpha_id, var]);
let exp_ax = arena.exp(alpha_x);
arena.mul(&[exp_ax, trig_part])
};
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
fn try_undetermined_trig_exp(
arena: &mut Arena,
rhs: ExprId,
b: &Q,
c: &Q,
var: ExprId,
) -> Option<ExprId> {
use num_traits::Zero;
let (coeff_r, term) = arena.as_coeff_term(rhs);
let node = arena.node(term).clone();
match node {
ExprNode::Sin(inner) => {
let (omega, constant) = extract_linear_numeric(arena, inner, var)?;
if !constant.is_zero() {
return None;
}
let zero_r = num_rational::Ratio::<num_bigint::BigInt>::zero();
try_trig_particular(arena, &coeff_r, &zero_r, &omega, b, c, var)
}
ExprNode::Cos(inner) => {
let (omega, constant) = extract_linear_numeric(arena, inner, var)?;
if !constant.is_zero() {
return None;
}
let zero_r = num_rational::Ratio::<num_bigint::BigInt>::zero();
try_trig_particular(arena, &zero_r, &coeff_r, &omega, b, c, var)
}
ExprNode::Exp(inner) => {
let (r_val, constant) = extract_linear_numeric(arena, inner, var)?;
if !constant.is_zero() {
return None;
}
try_exp_particular(arena, &coeff_r, &r_val, b, c, var)
}
_ => None,
}
}
fn extract_linear_numeric(arena: &Arena, expr: ExprId, var: ExprId) -> Option<(Q, Q)> {
let poly = crate::poly::polybridge::expr_to_poly(arena, expr, var)?;
if poly.degree()? != 1 {
return None;
}
Some((poly.coeff(1), poly.coeff(0)))
}
fn try_trig_particular(
arena: &mut Arena,
p: &Q,
q: &Q,
omega: &Q,
b: &Q,
c: &Q,
var: ExprId,
) -> Option<ExprId> {
use num_traits::Zero;
let omega_sq = omega * omega;
let d = c - &omega_sq; let bw = b * omega;
let det = &d * &d + &bw * &bw;
if !det.is_zero() {
let alpha = (&d * p + &bw * q) / &det;
let beta = (&d * q - &bw * p) / &det;
let omega_id = ode_ratio_to_expr(arena, omega);
let omega_x = arena.mul(&[omega_id, var]);
let mut terms = Vec::new();
if !alpha.is_zero() {
let alpha_id = ode_ratio_to_expr(arena, &alpha);
let sin_wx = arena.sin(omega_x);
terms.push(arena.mul(&[alpha_id, sin_wx]));
}
if !beta.is_zero() {
let beta_id = ode_ratio_to_expr(arena, &beta);
let cos_wx = arena.cos(omega_x);
terms.push(arena.mul(&[beta_id, cos_wx]));
}
match terms.len() {
0 => Some(arena.zero),
1 => Some(terms[0]),
_ => Some(arena.add(&terms)),
}
} else {
if omega.is_zero() {
return None;
}
let two_omega = num_rational::Ratio::from_integer(num_bigint::BigInt::from(2)) * omega;
let alpha = q / &two_omega;
let beta = -(p / &two_omega);
let omega_id = ode_ratio_to_expr(arena, omega);
let omega_x = arena.mul(&[omega_id, var]);
let mut inner_terms = Vec::new();
if !alpha.is_zero() {
let alpha_id = ode_ratio_to_expr(arena, &alpha);
let sin_wx = arena.sin(omega_x);
inner_terms.push(arena.mul(&[alpha_id, sin_wx]));
}
if !beta.is_zero() {
let beta_id = ode_ratio_to_expr(arena, &beta);
let cos_wx = arena.cos(omega_x);
inner_terms.push(arena.mul(&[beta_id, cos_wx]));
}
if inner_terms.is_empty() {
Some(arena.zero)
} else {
let inner = if inner_terms.len() == 1 {
inner_terms[0]
} else {
arena.add(&inner_terms)
};
Some(arena.mul(&[var, inner]))
}
}
}
fn try_exp_particular(
arena: &mut Arena,
coeff_r: &Q,
r: &Q,
b: &Q,
c: &Q,
var: ExprId,
) -> Option<ExprId> {
use num_traits::Zero;
let char_val = r * r + b * r + c;
let r_id = ode_ratio_to_expr(arena, r);
let rx = arena.mul(&[r_id, var]);
let exp_rx = arena.exp(rx);
if !char_val.is_zero() {
let a_val = coeff_r / &char_val;
let a_id = ode_ratio_to_expr(arena, &a_val);
Some(arena.mul(&[a_id, exp_rx]))
} else {
let deriv_val = num_rational::Ratio::from_integer(num_bigint::BigInt::from(2)) * r + b;
if !deriv_val.is_zero() {
let a_val = coeff_r / &deriv_val;
let a_id = ode_ratio_to_expr(arena, &a_val);
let x_exp = arena.mul(&[var, exp_rx]);
Some(arena.mul(&[a_id, x_exp]))
} else {
let two_r_val = num_rational::Ratio::from_integer(num_bigint::BigInt::from(2));
let a_val = coeff_r / &two_r_val;
let a_id = ode_ratio_to_expr(arena, &a_val);
let two_id = arena.int(2);
let x_sq = arena.pow(var, two_id);
let x2_exp = arena.mul(&[x_sq, exp_rx]);
Some(arena.mul(&[a_id, x2_exp]))
}
}
}
fn try_bernoulli(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying Bernoulli");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if expr_contains(arena, expr, d2y_dx2) {
return None;
}
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
use num_traits::Zero;
let mut dy_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut p_x_terms: Vec<ExprId> = Vec::new();
let mut q_x_terms: Vec<(ExprId, Q)> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
dy_coeff += coeff;
} else if !contains_sym(arena, child, func_sym) {
return None; } else if let Some(px) = extract_coeff_of_func(arena, child, func, func_sym, var_sym) {
p_x_terms.push(px);
} else if let Some((qx, n)) = extract_bernoulli_term(arena, child, func, func_sym, var_sym)
{
q_x_terms.push((qx, n));
} else {
return None;
}
}
if dy_coeff.is_zero() || q_x_terms.is_empty() {
return None;
}
let n_val = q_x_terms[0].1.clone();
if n_val.is_zero() || n_val.is_one() {
return None;
}
for (_, n) in &q_x_terms[1..] {
if *n != n_val {
return None;
}
}
let inv_dy = num_rational::Ratio::<num_bigint::BigInt>::one() / &dy_coeff;
let inv_dy_id = ode_ratio_to_expr(arena, &inv_dy);
let p_raw = if p_x_terms.is_empty() {
arena.zero
} else if p_x_terms.len() == 1 {
p_x_terms[0]
} else {
arena.add(&p_x_terms)
};
let p_x = if dy_coeff.is_one() {
p_raw
} else {
let s = arena.mul(&[inv_dy_id, p_raw]);
crate::transforms::eval::eval(arena, s)
};
let q_sum: Vec<ExprId> = q_x_terms.iter().map(|(qx, _)| *qx).collect();
let q_raw = if q_sum.len() == 1 {
q_sum[0]
} else {
arena.add(&q_sum)
};
let neg_q_raw = arena.neg(q_raw);
let q_x = if dy_coeff.is_one() {
neg_q_raw
} else {
let s = arena.mul(&[inv_dy_id, neg_q_raw]);
crate::transforms::eval::eval(arena, s)
};
if contains_sym(arena, p_x, func_sym) || contains_sym(arena, q_x, func_sym) {
return None;
}
let one_minus_n = num_rational::Ratio::<num_bigint::BigInt>::one() - &n_val;
let one_minus_n_id = ode_ratio_to_expr(arena, &one_minus_n);
let new_p = arena.mul(&[one_minus_n_id, p_x]);
let new_p = crate::transforms::eval::eval(arena, new_p);
let new_q = arena.mul(&[one_minus_n_id, q_x]);
let new_q = crate::transforms::eval::eval(arena, new_q);
let v = arena.symbol("__v");
let dv = arena.intern(ExprNode::Derivative(v, var));
let pv = arena.mul(&[new_p, v]);
let neg_new_q = arena.neg(new_q);
let linear_expr = arena.add(&[dv, pv, neg_new_q]);
let v_sym = match arena.node(v) {
ExprNode::Symbol(sid) => *sid,
_ => return None,
};
let v_result = try_first_order_linear_general(arena, linear_expr, v, var, v_sym, var_sym)
.or_else(|| try_simple_separable(arena, linear_expr, v, var, v_sym, var_sym))?;
let inv_one_minus_n = num_rational::Ratio::<num_bigint::BigInt>::one() / &one_minus_n;
let inv_id = ode_ratio_to_expr(arena, &inv_one_minus_n);
let solution = arena.pow(v_result.solution, inv_id);
let solution = crate::transforms::eval::eval(arena, solution);
let two = num_bigint::BigInt::from(2);
let branches = if inv_one_minus_n.denom() % &two == num_bigint::BigInt::from(0) {
let negative = arena.neg(solution);
vec![crate::transforms::eval::eval(arena, negative)]
} else if !inv_one_minus_n.is_integer() && inv_one_minus_n.numer().abs().is_one() {
let neg_v = arena.neg(v_result.solution);
let root = arena.pow(neg_v, inv_id);
let negative = arena.neg(root);
vec![crate::transforms::eval::eval(arena, negative)]
} else {
Vec::new()
};
Some(OdeResult {
solution,
constants: v_result.constants,
branches,
})
}
fn extract_bernoulli_term(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<(ExprId, Q)> {
if let ExprNode::Pow(base, exp) = arena.node(expr).clone()
&& base == func
{
let n = arena.as_num(exp)?.clone();
return Some((arena.one, n));
}
if let ExprNode::Mul(ref factors) = arena.node(expr).clone() {
let mut yn_idx = None;
let mut yn_exp = None;
for (i, &f) in factors.iter().enumerate() {
if let ExprNode::Pow(base, exp) = arena.node(f).clone()
&& base == func
&& let Some(n) = arena.as_num(exp)
{
yn_idx = Some(i);
yn_exp = Some(n.clone());
break;
}
}
if let (Some(idx), Some(n)) = (yn_idx, yn_exp) {
let mut other: Vec<ExprId> = Vec::new();
for (i, &f) in factors.iter().enumerate() {
if i != idx {
if contains_sym(arena, f, func_sym) {
return None;
}
other.push(f);
}
}
let qx = match other.len() {
0 => arena.one,
1 => other[0],
_ => arena.mul(&other),
};
return Some((qx, n));
}
}
{
let (coeff, term) = arena.as_coeff_term(expr);
if !coeff.is_one()
&& term != expr
&& let Some((inner_qx, n)) =
extract_bernoulli_term(arena, term, func, func_sym, _var_sym)
{
let coeff_id = ode_ratio_to_expr(arena, &coeff);
let qx = arena.mul(&[coeff_id, inner_qx]);
return Some((qx, n));
}
}
if let ExprNode::Neg(inner) = arena.node(expr).clone()
&& let Some((inner_qx, n)) = extract_bernoulli_term(arena, inner, func, func_sym, _var_sym)
{
let neg_qx = arena.neg(inner_qx);
return Some((neg_qx, n));
}
None
}
fn try_euler_cauchy(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying Euler-Cauchy");
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if !expr_contains(arena, expr, d2y_dx2) {
return None;
}
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
use num_traits::Zero;
let mut a_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut b_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut c_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let two_id = arena.int(2);
let x_sq = arena.pow(var, two_id);
let mut forcing: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == func {
c_coeff += coeff;
} else if !contains_sym(arena, child, func_sym) {
forcing.push(child);
} else if let ExprNode::Mul(ref factors) = arena.node(term).clone() {
let has_d2 = factors.contains(&d2y_dx2);
let has_d1 = factors.contains(&dy_dx);
let has_xsq = factors.contains(&x_sq);
let has_xvar = factors.contains(&var);
let other: Vec<ExprId> = factors
.iter()
.copied()
.filter(|&f| f != d2y_dx2 && f != dy_dx && f != x_sq && f != var)
.collect();
for &of in &other {
if contains_sym(arena, of, func_sym) || contains_sym(arena, of, var_sym) {
return None;
}
}
let extra = if other.is_empty() {
num_rational::Ratio::<num_bigint::BigInt>::from_integer(1.into())
} else {
let e = if other.len() == 1 {
other[0]
} else {
arena.mul(&other)
};
arena.as_num(e)?.clone()
};
if has_d2 && has_xsq && !has_d1 && !has_xvar {
a_coeff += &coeff * &extra;
} else if has_d1 && has_xvar && !has_d2 && !has_xsq {
b_coeff += &coeff * &extra;
} else {
return None;
}
} else {
return None;
}
}
if a_coeff.is_zero() {
return None;
}
if !forcing.is_empty() {
return euler_cauchy_forced(arena, &a_coeff, &b_coeff, &c_coeff, &forcing, var, var_sym);
}
let p = (&b_coeff - &a_coeff) / &a_coeff;
let q = &c_coeff / &a_coeff;
let four = num_rational::Ratio::<num_bigint::BigInt>::from_integer(4.into());
let disc = &p * &p - &four * &q;
let c1 = arena.symbol("C1");
let c2 = arena.symbol("C2");
if disc.is_positive() {
let r_var = arena.symbol("__r");
let r_sq = arena.pow(r_var, two_id);
let p_id = ode_ratio_to_expr(arena, &p);
let q_id = ode_ratio_to_expr(arena, &q);
let p_r = arena.mul(&[p_id, r_var]);
let char_eq = arena.add(&[r_sq, p_r, q_id]);
let roots = crate::transforms::solve::solve(arena, char_eq, r_var);
if roots.len() >= 2 {
let r1 = roots[0].value;
let r2 = roots[1].value;
let x_r1 = arena.pow(var, r1);
let x_r2 = arena.pow(var, r2);
let t1 = arena.mul(&[c1, x_r1]);
let t2 = arena.mul(&[c2, x_r2]);
let solution = arena.add(&[t1, t2]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
} else {
None
}
} else if disc.is_zero() {
let two_r = num_rational::Ratio::<num_bigint::BigInt>::from_integer(2.into());
let r = -(&p) / &two_r;
let r_id = ode_ratio_to_expr(arena, &r);
let x_r = arena.pow(var, r_id);
let ln_x = arena.ln(var);
let c2_ln = arena.mul(&[c2, ln_x]);
let inner = arena.add(&[c1, c2_ln]);
let solution = arena.mul(&[inner, x_r]);
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
} else {
let two_r = num_rational::Ratio::<num_bigint::BigInt>::from_integer(2.into());
let alpha = -(&p) / &two_r;
let neg_disc = -disc;
let neg_disc_id = ode_ratio_to_expr(arena, &neg_disc);
let half = arena.rational(1, 2);
let sqrt_neg_disc = arena.pow(neg_disc_id, half);
let two_expr = arena.int(2);
let beta = arena.div(sqrt_neg_disc, two_expr);
let beta = crate::transforms::eval::eval(arena, beta);
let ln_x = arena.ln(var);
let beta_ln_x = arena.mul(&[beta, ln_x]);
let cos_part = arena.cos(beta_ln_x);
let sin_part = arena.sin(beta_ln_x);
let c1_cos = arena.mul(&[c1, cos_part]);
let c2_sin = arena.mul(&[c2, sin_part]);
let trig_part = arena.add(&[c1_cos, c2_sin]);
let solution = if alpha.is_zero() {
trig_part
} else {
let alpha_id = ode_ratio_to_expr(arena, &alpha);
let x_alpha = arena.pow(var, alpha_id);
arena.mul(&[x_alpha, trig_part])
};
Some(OdeResult {
solution,
constants: vec![c1, c2],
branches: Vec::new(),
})
}
}
fn euler_cauchy_forced(
arena: &mut Arena,
a: &Q,
b: &Q,
c: &Q,
forcing: &[ExprId],
var: ExprId,
var_sym: SymbolId,
) -> Option<OdeResult> {
euler_cauchy_by_exp(arena, &[c.clone(), b - a, a.clone()], forcing, var, var_sym)
}
fn try_euler_cauchy_n(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
let chain = derivative_chain(arena, func, var);
let order = (1..chain.len())
.rev()
.find(|&k| expr_contains(arena, expr, chain[k]))?;
if order < 3 {
return None;
}
let ExprNode::Add(children) = arena.node(expr).clone() else {
return None;
};
let mut a = vec![Q::from_integer(0.into()); order + 1];
let mut forcing = Vec::new();
for &child in children.iter() {
if !contains_sym(arena, child, func_sym) {
forcing.push(child);
continue;
}
let (coeff, term) = arena.as_coeff_term(child);
let factors: Vec<ExprId> = match arena.node(term) {
ExprNode::Mul(fs) => fs.to_vec(),
_ => vec![term],
};
let mut k_deriv: Option<usize> = None;
let mut x_power: i64 = 0;
for f in factors {
if let Some(k) = chain.iter().position(|&d| d == f) {
if k_deriv.replace(k).is_some() {
return None;
}
} else if f == var {
x_power += 1;
} else if let ExprNode::Pow(b, e) = *arena.node(f)
&& b == var
&& let Some(p) = arena.as_num(e).filter(|p| p.is_integer())
{
x_power += num_traits::ToPrimitive::to_i64(&p.to_integer())?;
} else {
return None;
}
}
let k = k_deriv?;
if x_power != i64::try_from(k).ok()? {
return None;
}
a[k] += coeff;
}
if num_traits::Zero::is_zero(&a[order]) {
return None;
}
let mut cd = vec![Q::from_integer(0.into()); order + 1];
let mut falling = vec![Q::from_integer(1.into())];
for (k, ak) in a.iter().enumerate() {
if k > 0 {
let shift = Q::from_integer(((k - 1) as i64).into());
let mut next = vec![Q::from_integer(0.into()); falling.len() + 1];
for (j, v) in falling.iter().enumerate() {
next[j + 1] += v;
next[j] -= v * &shift;
}
falling = next;
}
for (j, v) in falling.iter().enumerate() {
cd[j] += ak * v;
}
}
euler_cauchy_by_exp(arena, &cd, &forcing, var, var_sym)
}
fn euler_cauchy_by_exp(
arena: &mut Arena,
cd: &[Q],
forcing: &[ExprId],
var: ExprId,
var_sym: SymbolId,
) -> Option<OdeResult> {
let order = cd.len().checked_sub(1)?;
let t = arena.symbol("__t");
let big_y = arena.symbol("__Y");
let (t_sym, y_sym) = match (arena.node(t), arena.node(big_y)) {
(ExprNode::Symbol(ts), ExprNode::Symbol(ys)) => (*ts, *ys),
_ => return None,
};
let exp_t = arena.exp(t);
let ln_x = arena.ln(var);
let mut parts = Vec::with_capacity(cd.len() + forcing.len());
let mut d = big_y;
for (j, c) in cd.iter().enumerate() {
if j > 0 {
d = arena.intern(ExprNode::Derivative(d, t));
}
if num_traits::Zero::is_zero(c) {
continue;
}
let c_id = ode_ratio_to_expr(arena, c);
parts.push(arena.mul(&[c_id, d]));
}
for &f in forcing {
let g = crate::transforms::subs::subs(arena, f, ln_x, t);
let g = crate::transforms::subs::subs(arena, g, var, exp_t);
parts.push(crate::transforms::eval::eval(arena, g));
}
let ode_t = arena.add(&parts);
let ode_t = crate::transforms::eval::eval(arena, ode_t);
if contains_sym(arena, ode_t, var_sym) {
return None;
}
let res = dsolve(arena, ode_t, big_y, t)?;
if contains_sym(arena, res.solution, y_sym)
|| res.constants.len() != order
|| crate::base::walk::has_unevaluated(arena, res.solution)
{
return None;
}
let mut replacements: Vec<(ExprId, ExprId)> = Vec::new();
for id in crate::base::walk::post_order_ids(arena, res.solution) {
let ExprNode::Exp(arg) = arena.node(id).clone() else {
continue;
};
let arg = crate::transforms::expand::expand(arena, arg);
let arg = crate::transforms::eval::eval(arena, arg);
let Some(coeffs) = arena.coefficients_of(arg, t) else {
continue;
};
if coeffs.len() != 2 || contains_sym(arena, coeffs[1], t_sym) {
continue;
}
let x_k = arena.pow(var, coeffs[1]);
let rest = arena.exp(coeffs[0]);
let r = arena.mul(&[x_k, rest]);
let r = crate::transforms::eval::eval(arena, r);
replacements.push((id, r));
}
let solution = crate::transforms::subs::subs_map(arena, res.solution, &replacements);
let solution = crate::transforms::subs::subs(arena, solution, t, ln_x);
let solution = crate::transforms::eval::eval(arena, solution);
if contains_sym(arena, solution, t_sym) {
return None;
}
Some(OdeResult {
solution,
constants: res.constants,
branches: Vec::new(),
})
}
fn try_variation_of_parameters(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
tracing::debug!("ode: trying variation of parameters");
use num_traits::Zero;
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
let children = match arena.node(expr).clone() {
ExprNode::Add(children) => children,
_ => return None,
};
let mut a_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut b_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut c_coeff = num_rational::Ratio::<num_bigint::BigInt>::zero();
let mut f_terms: Vec<ExprId> = Vec::new();
for &child in &children {
let (coeff, term) = arena.as_coeff_term(child);
if term == d2y_dx2 {
a_coeff += coeff;
} else if term == dy_dx {
b_coeff += coeff;
} else if term == func {
c_coeff += coeff;
} else if !contains_sym(arena, child, func_sym) {
f_terms.push(child);
} else {
return None;
}
}
if a_coeff.is_zero() || f_terms.is_empty() {
return None;
}
let b = &b_coeff / &a_coeff;
let c = &c_coeff / &a_coeff;
let f_sum = if f_terms.len() == 1 {
f_terms[0]
} else {
arena.add(&f_terms)
};
let neg_f = arena.neg(f_sum);
let a_id = ode_ratio_to_expr(arena, &a_coeff);
let g_x = arena.div(neg_f, a_id);
let g_x = crate::transforms::eval::eval(arena, g_x);
let homo = solve_characteristic_equation(arena, b.clone(), c.clone(), var)?;
let (y1, y2) = extract_fundamental_solutions(arena, homo.solution, &homo.constants)?;
let y1_prime = crate::transforms::diff::diff(arena, y1, var);
let y2_prime = crate::transforms::diff::diff(arena, y2, var);
let w_term1 = arena.mul(&[y1, y2_prime]);
let w_term2 = arena.mul(&[y2, y1_prime]);
let wronskian = arena.sub(w_term1, w_term2);
let wronskian = crate::transforms::eval::eval(arena, wronskian);
let wronskian = crate::transforms::expand::expand(arena, wronskian);
let wronskian = crate::transforms::eval::eval(arena, wronskian);
let wronskian = if contains_sym(arena, wronskian, var_sym) {
let w0 = crate::transforms::subs::subs(arena, wronskian, var, arena.zero);
let w0 = crate::transforms::eval::eval(arena, w0);
if w0 == arena.zero {
return None;
}
if b.is_zero() {
w0
} else {
let neg_b_id = ode_ratio_to_expr(arena, &(-b.clone()));
let neg_bx = arena.mul(&[neg_b_id, var]);
let exp_nbx = arena.exp(neg_bx);
arena.mul(&[w0, exp_nbx])
}
} else {
wronskian
};
if wronskian == arena.zero {
return None;
}
let y2_g = arena.mul(&[y2, g_x]);
let integrand1 = arena.div(y2_g, wronskian);
let integrand1 = crate::transforms::eval::eval(arena, integrand1);
let integral1 = integrate_forms(arena, integrand1, var);
if integral_failed(arena, integral1) {
return None;
}
let y1_g = arena.mul(&[y1, g_x]);
let integrand2 = arena.div(y1_g, wronskian);
let integrand2 = crate::transforms::eval::eval(arena, integrand2);
let integral2 = integrate_forms(arena, integrand2, var);
if integral_failed(arena, integral2) {
return None;
}
let term1 = arena.mul(&[y1, integral1]);
let neg_term1 = arena.neg(term1);
let term2 = arena.mul(&[y2, integral2]);
let y_p = arena.add(&[neg_term1, term2]);
let y_p = crate::transforms::eval::eval(arena, y_p);
let solution = arena.add(&[homo.solution, y_p]);
let solution = crate::transforms::eval::eval(arena, solution);
Some(OdeResult {
solution,
constants: homo.constants,
branches: Vec::new(),
})
}
fn extract_fundamental_solutions(
arena: &mut Arena,
homo_solution: ExprId,
constants: &[ExprId],
) -> Option<(ExprId, ExprId)> {
if constants.len() != 2 {
return None;
}
let c1 = constants[0];
let c2 = constants[1];
let zero = arena.zero;
let one = arena.one;
let y1 = crate::transforms::subs::subs(arena, homo_solution, c1, one);
let y1 = crate::transforms::subs::subs(arena, y1, c2, zero);
let y1 = crate::transforms::eval::eval(arena, y1);
let y2 = crate::transforms::subs::subs(arena, homo_solution, c1, zero);
let y2 = crate::transforms::subs::subs(arena, y2, c2, one);
let y2 = crate::transforms::eval::eval(arena, y2);
if y1 == arena.zero || y2 == arena.zero {
return None;
}
Some((y1, y2))
}
const MAX_ODE_ORDER: usize = 12;
struct LinearCcOde {
coeffs: Vec<Q>,
forcing: Vec<ExprId>,
}
fn derivative_chain(arena: &mut Arena, func: ExprId, var: ExprId) -> Vec<ExprId> {
let mut chain = vec![func];
for k in 1..=MAX_ODE_ORDER {
let d = arena.intern(ExprNode::Derivative(chain[k - 1], var));
chain.push(d);
}
chain
}
fn extract_linear_cc(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<LinearCcOde> {
use num_traits::Zero;
let chain = derivative_chain(arena, func, var);
let children: Vec<ExprId> = match arena.node(expr).clone() {
ExprNode::Add(c) => c.to_vec(),
_ => vec![expr],
};
let mut coeffs = vec![Q::zero(); MAX_ODE_ORDER + 1];
let mut forcing = Vec::new();
for child in children {
let (coeff, term) = arena.as_coeff_term(child);
if let Some(k) = chain.iter().position(|&d| d == term) {
coeffs[k] += coeff;
} else if !contains_sym(arena, child, func_sym) {
forcing.push(child);
} else {
return None;
}
}
while coeffs.len() > 1 && coeffs.last().is_some_and(|c| c.is_zero()) {
coeffs.pop();
}
if coeffs.len() < 2 {
return None; }
Some(LinearCcOde { coeffs, forcing })
}
#[derive(Clone, Debug, PartialEq)]
struct CQ {
re: Q,
im: Q,
}
impl CQ {
fn new(re: Q, im: Q) -> Self {
Self { re, im }
}
fn real(re: Q) -> Self {
Self {
re,
im: num_traits::Zero::zero(),
}
}
fn is_zero(&self) -> bool {
num_traits::Zero::is_zero(&self.re) && num_traits::Zero::is_zero(&self.im)
}
fn add(&self, o: &CQ) -> CQ {
CQ::new(&self.re + &o.re, &self.im + &o.im)
}
fn sub(&self, o: &CQ) -> CQ {
CQ::new(&self.re - &o.re, &self.im - &o.im)
}
fn mul(&self, o: &CQ) -> CQ {
CQ::new(
&self.re * &o.re - &self.im * &o.im,
&self.re * &o.im + &self.im * &o.re,
)
}
fn div(&self, o: &CQ) -> CQ {
let denom = &o.re * &o.re + &o.im * &o.im;
let num = self.mul(&CQ::new(o.re.clone(), -o.im.clone()));
CQ::new(num.re / &denom, num.im / &denom)
}
fn scale(&self, r: &Q) -> CQ {
CQ::new(&self.re * r, &self.im * r)
}
}
fn shifted_coefficient(p: &[Q], j: usize, lambda: &CQ) -> CQ {
let mut acc = CQ::real(num_traits::Zero::zero());
let mut lambda_pow = CQ::real(num_traits::One::one());
for (k, a_k) in p.iter().enumerate().skip(j) {
let binom = Q::from_integer(binomial(k as u64, j as u64));
let term = lambda_pow.scale(&(a_k * binom));
acc = acc.add(&term);
lambda_pow = lambda_pow.mul(lambda);
}
acc
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
enum TrigKind {
None,
Cos,
Sin,
}
struct ForcingTerm {
coeff: Q,
degree: usize,
a: Q,
b: Q,
kind: TrigKind,
}
fn parse_forcing_term(arena: &mut Arena, term: ExprId, var: ExprId) -> Option<ForcingTerm> {
use num_traits::{Signed, Zero};
let (mut coeff, t) = arena.as_coeff_term(term);
let factors: Vec<ExprId> = if t == arena.one {
Vec::new()
} else {
match arena.node(t).clone() {
ExprNode::Mul(c) => c.to_vec(),
_ => vec![t],
}
};
let mut degree = 0usize;
let mut a = Q::zero();
let mut b = Q::zero();
let mut kind = TrigKind::None;
let mut seen_exp = false;
for f in factors {
if f == var {
degree += 1;
continue;
}
match arena.node(f).clone() {
ExprNode::Pow(base, exp) if base == var => {
let n = arena.as_num(exp)?.clone();
if !n.is_integer() || n.is_negative() {
return None;
}
degree += usize::try_from(n.to_integer()).ok()?;
}
ExprNode::Exp(inner) => {
if seen_exp {
return None;
}
let (rate, c) = extract_linear_numeric(arena, inner, var)?;
if !c.is_zero() {
return None;
}
seen_exp = true;
a = rate;
}
ExprNode::Sin(inner) => {
if kind != TrigKind::None {
return None;
}
let (freq, c) = extract_linear_numeric(arena, inner, var)?;
if !c.is_zero() || freq.is_zero() {
return None;
}
kind = TrigKind::Sin;
if freq.is_negative() {
coeff = -coeff;
b = -freq;
} else {
b = freq;
}
}
ExprNode::Cos(inner) => {
if kind != TrigKind::None {
return None;
}
let (freq, c) = extract_linear_numeric(arena, inner, var)?;
if !c.is_zero() || freq.is_zero() {
return None;
}
kind = TrigKind::Cos;
b = freq.abs();
}
_ => {
let r = arena.as_num(f)?;
coeff *= r.clone();
}
}
}
Some(ForcingTerm {
coeff,
degree,
a,
b,
kind,
})
}
fn rat_poly_expr(arena: &mut Arena, coeffs: &[Q], var: ExprId) -> ExprId {
build_polynomial_expr(arena, coeffs, var)
}
fn particular_for_group(
arena: &mut Arena,
p: &[Q],
q: &[Q],
a: &Q,
b: &Q,
kind: TrigKind,
var: ExprId,
) -> Option<ExprId> {
use num_traits::Zero;
let lambda = CQ::new(a.clone(), b.clone());
let n = p.len() - 1;
let c: Vec<CQ> = (0..=n)
.map(|j| shifted_coefficient(p, j, &lambda))
.collect();
let s = c.iter().position(|cj| !cj.is_zero())?;
let d = q.len().checked_sub(1)?;
let m = d + s;
let mut w = vec![CQ::real(Q::zero()); m + 1];
for k in (0..=d).rev() {
let mut rhs = CQ::real(q[k].clone());
for j in (s + 1)..=n {
if k + j > m {
continue;
}
let fact = falling_factorial_rat(k + j, j);
rhs = rhs.sub(&c[j].mul(&w[k + j]).scale(&fact));
}
let lead = c[s].scale(&falling_factorial_rat(k + s, s));
w[k + s] = rhs.div(&lead);
}
let wr: Vec<Q> = w.iter().map(|z| z.re.clone()).collect();
let wi: Vec<Q> = w.iter().map(|z| z.im.clone()).collect();
let wr_x = rat_poly_expr(arena, &wr, var);
let wi_x = rat_poly_expr(arena, &wi, var);
let exp_ax = if a.is_zero() {
arena.one
} else {
let a_id = ode_ratio_to_expr(arena, a);
let ax = arena.mul(&[a_id, var]);
arena.exp(ax)
};
let body = match kind {
TrigKind::None => wr_x,
TrigKind::Cos | TrigKind::Sin => {
let b_id = ode_ratio_to_expr(arena, b);
let bx = arena.mul(&[b_id, var]);
let cos_bx = arena.cos(bx);
let sin_bx = arena.sin(bx);
if kind == TrigKind::Cos {
let t1 = arena.mul(&[wr_x, cos_bx]);
let t2 = arena.mul(&[wi_x, sin_bx]);
arena.sub(t1, t2)
} else {
let t1 = arena.mul(&[wr_x, sin_bx]);
let t2 = arena.mul(&[wi_x, cos_bx]);
arena.add(&[t1, t2])
}
}
};
let y_p = arena.mul(&[exp_ax, body]);
let y_p = crate::transforms::expand::expand(arena, y_p);
Some(crate::transforms::eval::eval(arena, y_p))
}
fn falling_factorial_rat(n: usize, j: usize) -> Q {
let mut r = Q::from_integer(1.into());
for i in 0..j {
r *= Q::from_integer(((n - i) as i64).into());
}
r
}
fn homogeneous_basis_cc(arena: &mut Arena, coeffs: &[Q], var: ExprId) -> Option<Vec<ExprId>> {
let p = crate::poly::Poly::from_coeffs(coeffs.to_vec());
let r_sym = arena.symbol("__r_cc");
let mut basis: Vec<ExprId> = Vec::new();
for (factor, mult) in p.squarefree_factors() {
for root in characteristic_roots(arena, &factor, r_sym, var)? {
for j in 0..mult {
match root {
CharRoot::Real(r) => basis.push(mode_real(arena, r, j, var)),
CharRoot::Pair { re, im } => {
let (c, s) = mode_complex(arena, re, im, j, var);
basis.push(c);
basis.push(s);
}
}
}
}
}
Some(basis)
}
#[derive(Clone, Copy, Debug)]
enum CharRoot {
Real(ExprId),
Pair { re: ExprId, im: ExprId },
}
fn characteristic_roots(
arena: &mut Arena,
f: &crate::poly::Poly,
r_sym: ExprId,
var: ExprId,
) -> Option<Vec<CharRoot>> {
if let Some(roots) = radical_char_roots(arena, f, r_sym) {
return Some(roots);
}
let (_, irreducible) = f.factor_over_z();
let mut out = Vec::new();
for (g, _) in &irreducible {
let radical = if irreducible.len() > 1 {
radical_char_roots(arena, g, r_sym)
} else {
None
};
let roots = match radical {
Some(r) => r,
None => rootof_char_roots(arena, g, var)?,
};
out.extend(roots);
}
Some(out)
}
const IM_REAL: f64 = 1e-9;
const IM_NONREAL: f64 = 1e-6;
fn radical_char_roots(
arena: &mut Arena,
f: &crate::poly::Poly,
r_sym: ExprId,
) -> Option<Vec<CharRoot>> {
let deg = f.degree()?;
let f_expr = crate::poly::polybridge::poly_to_expr(arena, f, r_sym);
let sols: Vec<ExprId> = crate::transforms::solve::solve(arena, f_expr, r_sym)
.into_iter()
.map(|s| s.value)
.collect();
if sols.len() != deg {
return None;
}
let mut values = Vec::with_capacity(deg);
for &s in &sols {
if crate::base::walk::contains(arena, s, r_sym) || has_rootof(arena, s) {
return None;
}
let z = crate::transforms::evalf::evalf_complex64(arena, s).ok()?;
if !(z.re.is_finite() && z.im.is_finite()) {
return None;
}
values.push(z);
}
let n_real = crate::poly::sturm::SturmChain::new(f).count_real_roots();
let rel_im = |z: &num_complex::Complex64| z.im.abs() / z.norm().max(1.0);
let mut order: Vec<usize> = (0..deg).collect();
order.sort_by(|&a, &b| rel_im(&values[a]).total_cmp(&rel_im(&values[b])));
if (n_real > 0 && rel_im(&values[order[n_real - 1]]) > IM_REAL)
|| (n_real < deg && rel_im(&values[order[n_real]]) < IM_NONREAL)
{
return None;
}
let mut is_real = vec![false; deg];
for &k in &order[..n_real] {
is_real[k] = true;
}
let close = |a: f64, b: f64| (a - b).abs() <= 1e-9 * (1.0 + b.abs());
let i_unit = arena.i_unit;
let mut used = vec![false; deg];
let mut out = Vec::new();
for i in 0..deg {
if used[i] {
continue;
}
used[i] = true;
if is_real[i] {
let mut r = sols[i];
if crate::base::walk::contains(arena, r, i_unit) {
r = real_imag_parts(arena, r).0;
if !real_radical_form(arena, r) {
return None;
}
let z = crate::transforms::evalf::evalf_complex64(arena, r).ok()?;
if !close(z.re, values[i].re) || z.im.abs() > IM_REAL * (1.0 + z.re.abs()) {
return None;
}
}
out.push(CharRoot::Real(r));
continue;
}
let conj = values[i].conj();
let j = (0..deg)
.filter(|&k| !used[k] && !is_real[k])
.min_by(|&a, &b| {
(values[a] - conj)
.norm()
.total_cmp(&(values[b] - conj).norm())
})?;
if (values[j] - conj).norm() > 1e-9 * (1.0 + conj.norm()) {
return None;
}
used[j] = true;
let pos = if values[i].im > 0.0 { i } else { j };
let (re, im) = real_imag_parts(arena, sols[pos]);
if !real_radical_form(arena, re) || !real_radical_form(arena, im) {
return None;
}
let zr = crate::transforms::evalf::evalf_complex64(arena, re).ok()?;
let zi = crate::transforms::evalf::evalf_complex64(arena, im).ok()?;
if !close(zr.re, values[pos].re) || !close(zi.re, values[pos].im) {
return None;
}
out.push(CharRoot::Pair { re, im });
}
Some(out)
}
fn real_imag_parts(arena: &mut Arena, z: ExprId) -> (ExprId, ExprId) {
let (re, im) = arena.as_real_imag_expr(z);
let mut parts = [re, im];
for part in &mut parts {
let post = crate::base::walk::post_order_ids(arena, *part);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match arena.node(rebuilt).clone() {
ExprNode::Re(u) if is_real_radical(arena, u) => u,
ExprNode::Im(u) if is_real_radical(arena, u) => arena.zero,
ExprNode::Cos(a) | ExprNode::Sin(a) => {
let (c, t) = arena.as_coeff_term(a);
if c.is_negative() && t == arena.pi {
let pos = ode_ratio_to_expr(arena, &(-c));
let pos = arena.mul(&[pos, t]);
if matches!(arena.node(rebuilt), ExprNode::Cos(_)) {
arena.cos(pos)
} else {
let s = arena.sin(pos);
arena.neg(s)
}
} else {
rebuilt
}
}
_ => rebuilt,
};
cache.insert(id, new);
}
let out = cache.get(part).copied().unwrap_or(*part);
*part = crate::transforms::eval::eval(arena, out);
}
(parts[0], parts[1])
}
fn is_real_radical(arena: &Arena, u: ExprId) -> bool {
if !real_radical_form(arena, u) || !crate::base::walk::free_symbols(arena, u).is_empty() {
return false;
}
let positive = |b: ExprId| {
crate::transforms::evalf::evalf_complex64(arena, b)
.is_ok_and(|z| z.im == 0.0 && z.re > 1e-12 * (1.0 + z.re.abs()))
};
let mut stack = vec![u];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
let node = arena.node(n);
match node {
ExprNode::Pow(b, e) if !arena.as_num(*e).is_some_and(Q::is_integer) => {
if !positive(*b) {
return false;
}
}
ExprNode::Ln(b) if !positive(*b) => return false,
_ => {}
}
node.for_each_child(|c| stack.push(c));
}
crate::transforms::evalf::evalf_complex64(arena, u)
.is_ok_and(|z| z.im == 0.0 && z.re.is_finite())
}
fn real_radical_form(arena: &Arena, e: ExprId) -> bool {
let i_unit = arena.i_unit;
let mut stack = vec![e];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
if n == i_unit {
return false;
}
let node = arena.node(n);
if matches!(
node,
ExprNode::Re(_) | ExprNode::Im(_) | ExprNode::Arg(_) | ExprNode::Conjugate(_)
) {
return false;
}
node.for_each_child(|c| stack.push(c));
}
true
}
fn horner_complex(
arena: &mut Arena,
h: &crate::poly::Poly,
re: ExprId,
im: ExprId,
) -> (ExprId, ExprId) {
let (mut a, mut b) = (arena.zero, arena.zero);
for c in h.coeffs().iter().rev() {
let a_re = arena.mul(&[a, re]);
let b_im = arena.mul(&[b, im]);
let neg_b_im = arena.neg(b_im);
let c = ode_ratio_to_expr(arena, c);
let a_new = arena.add(&[a_re, neg_b_im, c]);
let a_im = arena.mul(&[a, im]);
let b_re = arena.mul(&[b, re]);
let b_new = arena.add(&[a_im, b_re]);
a = crate::transforms::eval::eval(arena, a_new);
b = crate::transforms::eval::eval(arena, b_new);
}
(a, b)
}
fn has_rootof(arena: &Arena, e: ExprId) -> bool {
let mut stack = vec![e];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(n) = stack.pop() {
if !seen.insert(n) {
continue;
}
let node = arena.node(n);
if matches!(node, ExprNode::RootOf(..)) {
return true;
}
node.for_each_child(|c| stack.push(c));
}
false
}
const ROOTOF_PREC: usize = 128;
fn rootof_char_roots(
arena: &mut Arena,
g: &crate::poly::Poly,
var: ExprId,
) -> Option<Vec<CharRoot>> {
let deg = g.degree()?;
let balls = crate::poly::roots::certified_roots(g, ROOTOF_PREC);
if balls.len() != deg
|| balls.iter().any(|b| b.radius.is_none())
|| balls.iter().filter(|b| b.real).count()
!= crate::poly::sturm::SturmChain::new(g).count_real_roots()
{
return None;
}
let var_name = match arena.node(var) {
ExprNode::Symbol(s) => arena.symbol_name(*s).to_string(),
_ => String::new(),
};
let name = ["r", "s", "u", "w"]
.into_iter()
.find(|n| *n != var_name)
.unwrap_or("r");
let bound = arena.symbol(name);
let g_expr = crate::poly::polybridge::poly_to_expr(arena, g, bound);
let mut out = Vec::new();
let mut pairs = 0usize;
for (k, ball) in balls.iter().enumerate() {
let idx = arena.int(i64::try_from(k).ok()?);
let root = arena.intern(ExprNode::RootOf(g_expr, bound, idx));
if ball.real {
out.push(CharRoot::Real(root));
} else if ball.value.1.is_positive() {
let re = arena.re(root);
let im = arena.im(root);
out.push(CharRoot::Pair { re, im });
pairs += 1;
}
}
(out.len() + pairs == deg).then_some(out)
}
pub(crate) fn is_linear_autonomous(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
) -> bool {
let ExprNode::Symbol(func_sym) = *arena.node(func) else {
return false;
};
let chain = derivative_chain(arena, func, var);
let children: Vec<ExprId> = match arena.node(expr) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![expr],
};
let mut has_func = false;
for child in children {
let (_, term) = arena.as_coeff_term(child);
if symbolic_linear_term(arena, term, &chain, var, func_sym).is_some() {
has_func = true;
} else if contains_sym(arena, child, func_sym) {
return false;
}
}
has_func
}
pub(crate) fn is_cauchy_euler(arena: &mut Arena, expr: ExprId, func: ExprId, var: ExprId) -> bool {
let ExprNode::Symbol(func_sym) = *arena.node(func) else {
return false;
};
let ExprNode::Symbol(var_sym) = *arena.node(var) else {
return false;
};
let chain = derivative_chain(arena, func, var);
let children: Vec<ExprId> = match arena.node(expr) {
ExprNode::Add(c) => c.to_vec(),
_ => vec![expr],
};
let mut has_func = false;
for child in children {
if !contains_sym(arena, child, func_sym) {
continue;
}
let factors: Vec<ExprId> = match arena.node(child) {
ExprNode::Mul(fs) => fs.to_vec(),
ExprNode::Neg(inner) => vec![arena.neg_one, *inner],
_ => vec![child],
};
let mut order = None;
let mut rest = Vec::with_capacity(factors.len());
for f in factors {
match chain.iter().position(|&d| d == f) {
Some(k) if order.is_none() => order = Some(k),
Some(_) => return false,
None if contains_sym(arena, f, func_sym) => return false,
None => rest.push(f),
}
}
let Some(k) = order else {
return false;
};
let r = arena.mul(&rest);
let xk = {
let e = arena.int(i64::try_from(k).unwrap_or(0));
arena.pow(var, e)
};
let q = arena.div(r, xk);
let q = crate::transforms::eval::eval(arena, q);
if contains_sym(arena, q, var_sym) {
return false;
}
has_func = true;
}
has_func
}
pub(crate) fn linear_cc_ivp_lagrange(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
x0: ExprId,
values: &[ExprId],
) -> Option<ExprId> {
use crate::poly::Poly;
let ExprNode::Symbol(func_sym) = *arena.node(func) else {
return None;
};
let ode = extract_linear_cc(arena, expr, func, var, func_sym)?;
let n = ode.coeffs.len() - 1;
if n < 2 || values.len() != n {
return None;
}
let p = Poly::from_coeffs(ode.coeffs.clone()).make_monic();
if p.squarefree_factors().iter().any(|(_, m)| *m != 1) {
return None;
}
let r_sym = arena.symbol("__r_cc");
let (_, irreducible) = p.factor_over_z();
if irreducible.iter().all(|(g, _)| g.degree().unwrap_or(0) < 3) {
return None;
}
let mut groups: Vec<(Poly, Vec<CharRoot>)> = Vec::with_capacity(irreducible.len());
for (g, _) in irreducible {
let roots = match radical_char_roots(arena, &g, r_sym) {
Some(r) => r,
None => rootof_char_roots(arena, &g, var)?,
};
groups.push((g, roots));
}
let general = try_nth_order_linear_const_coeff(arena, expr, func, var, func_sym)?;
let mut y_p = general.solution;
let zero = arena.zero;
for &c in &general.constants {
y_p = crate::transforms::subs::subs(arena, y_p, c, zero);
}
let y_p = crate::transforms::eval::eval(arena, y_p);
let mut w = Vec::with_capacity(n);
let mut d = y_p;
for &v in values {
let at = crate::transforms::subs::subs(arena, d, var, x0);
let wk = arena.sub(v, at);
w.push(crate::transforms::eval::eval(arena, wk));
d = crate::transforms::diff::diff(arena, d, var);
}
let t_poly = Poly::from_coeffs(vec![Q::from_integer(0.into()), Q::from_integer(1.into())]);
let mut b = vec![Poly::from_coeffs(vec![Q::from_integer(1.into())]); n];
for k in (1..n).rev() {
let ak = Poly::from_coeffs(vec![p.coeff(k)]);
b[k - 1] = ak.add(&t_poly.mul(&b[k]));
}
let dp = p.derivative();
let t_sym = arena.symbol("__t_cc");
let shift = arena.sub(var, x0);
let two = arena.int(2);
let mut terms = vec![y_p];
for (g, roots) in &groups {
let e = Poly::extended_gcd(&dp.rem(g), g);
if e.gcd.degree() != Some(0) {
return None;
}
let inv_dp = e.x;
let h: Vec<Poly> = b.iter().map(|bk| bk.mul(&inv_dp).rem(g)).collect();
for root in roots {
match *root {
CharRoot::Real(r) => {
let mut c = Vec::with_capacity(n);
for (hk, wk) in h.iter().zip(&w) {
let hk = crate::poly::polybridge::poly_to_expr(arena, hk, t_sym);
let at = crate::transforms::subs::subs(arena, hk, t_sym, r);
c.push(arena.mul(&[*wk, at]));
}
let c = arena.add(&c);
let rx = arena.mul(&[r, shift]);
let e = arena.exp(rx);
terms.push(arena.mul(&[c, e]));
}
CharRoot::Pair { re, im } => {
let (mut a, mut bb) = (Vec::with_capacity(n), Vec::with_capacity(n));
for (hk, wk) in h.iter().zip(&w) {
let (h_re, h_im) = horner_complex(arena, hk, re, im);
a.push(arena.mul(&[*wk, h_re]));
bb.push(arena.mul(&[*wk, h_im]));
}
let a = arena.add(&a);
let a = arena.mul(&[two, a]);
let bsum = arena.add(&bb);
let neg_two = arena.int(-2);
let bcoef = arena.mul(&[neg_two, bsum]);
let ax = arena.mul(&[re, shift]);
let env = arena.exp(ax);
let bx = arena.mul(&[im, shift]);
let cos_bx = arena.cos(bx);
let sin_bx = arena.sin(bx);
let ca = arena.mul(&[a, cos_bx]);
let sb = arena.mul(&[bcoef, sin_bx]);
let osc = arena.add(&[ca, sb]);
terms.push(arena.mul(&[env, osc]));
}
}
}
}
let sol = arena.add(&terms);
Some(crate::transforms::eval::eval(arena, sol))
}
fn mode_real(arena: &mut Arena, r: ExprId, j: usize, var: ExprId) -> ExprId {
let x_pow = match j {
0 => arena.one,
1 => var,
_ => {
let e = arena.int(j as i64);
arena.pow(var, e)
}
};
if arena.is_zero_structural(r) {
return x_pow;
}
let rx = arena.mul(&[r, var]);
let exp_rx = arena.exp(rx);
let m = arena.mul(&[x_pow, exp_rx]);
crate::transforms::eval::eval(arena, m)
}
fn mode_complex(
arena: &mut Arena,
alpha: ExprId,
beta: ExprId,
j: usize,
var: ExprId,
) -> (ExprId, ExprId) {
let envelope = mode_real(arena, alpha, j, var);
let bx = arena.mul(&[beta, var]);
let cos_bx = arena.cos(bx);
let sin_bx = arena.sin(bx);
let c = arena.mul(&[envelope, cos_bx]);
let s = arena.mul(&[envelope, sin_bx]);
(
crate::transforms::eval::eval(arena, c),
crate::transforms::eval::eval(arena, s),
)
}
fn try_nth_order_linear_const_coeff(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<OdeResult> {
use num_traits::Zero;
let ode = extract_linear_cc(arena, expr, func, var, func_sym)?;
let n = ode.coeffs.len() - 1;
if n < 2 {
return None;
}
tracing::debug!(order = n, "ode: nth-order linear constant-coefficient");
let mut groups: Vec<((Q, Q, TrigKind), Vec<Q>)> = Vec::new();
for &term in &ode.forcing {
let ft = parse_forcing_term(arena, term, var)?;
let key = (ft.a.clone(), ft.b.clone(), ft.kind);
let entry = match groups.iter_mut().find(|(k, _)| *k == key) {
Some(e) => e,
None => {
groups.push((key, Vec::new()));
groups.last_mut()?
}
};
if entry.1.len() <= ft.degree {
entry.1.resize(ft.degree + 1, Q::zero());
}
entry.1[ft.degree] -= ft.coeff; }
let basis = homogeneous_basis_cc(arena, &ode.coeffs, var)?;
if basis.len() != n {
return None;
}
let mut constants = Vec::with_capacity(n);
let mut terms = Vec::with_capacity(n + groups.len());
for (i, &b) in basis.iter().enumerate() {
let c = arena.symbol(&format!("C{}", i + 1));
constants.push(c);
terms.push(arena.mul(&[c, b]));
}
for ((a, b, kind), q) in &groups {
if q.iter().all(Zero::is_zero) {
continue;
}
let y_p = particular_for_group(arena, &ode.coeffs, q, a, b, *kind, var)?;
terms.push(y_p);
}
let solution = arena.add(&terms);
let solution = crate::transforms::eval::eval(arena, solution);
Some(OdeResult {
solution,
constants,
branches: Vec::new(),
})
}
fn symbolic_linear_term(
arena: &mut Arena,
term: ExprId,
chain: &[ExprId],
var: ExprId,
func_sym: SymbolId,
) -> Option<(usize, ExprId)> {
if let Some(k) = chain.iter().position(|&d| d == term) {
return Some((k, arena.one));
}
let ExprNode::Mul(factors) = arena.node(term).clone() else {
return None;
};
let mut order = None;
let mut rest = Vec::with_capacity(factors.len());
for f in factors {
if let Some(k) = chain.iter().position(|&d| d == f) {
if order.replace(k).is_some() {
return None;
}
} else if crate::base::walk::contains(arena, f, var) || contains_sym(arena, f, func_sym) {
return None;
} else {
rest.push(f);
}
}
Some((order?, arena.mul(&rest)))
}
fn generic_branches(arena: &mut Arena, e: ExprId, var: ExprId) -> Option<ExprId> {
fn is_generic(arena: &Arena, c: ExprId) -> bool {
match arena.node(c) {
ExprNode::Ne(..) => true,
ExprNode::And(cs) => cs
.iter()
.all(|&k| matches!(arena.node(k), ExprNode::Ne(..))),
_ => false,
}
}
let post = crate::base::walk::post_order_ids(arena, e);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match arena.node(rebuilt).clone() {
ExprNode::Piecewise(pairs) => {
let (v, c) = *pairs.first()?;
if crate::base::walk::contains(arena, c, var) || !is_generic(arena, c) {
return None;
}
v
}
_ => rebuilt,
};
cache.insert(id, new);
}
let r = cache.get(&e).copied().unwrap_or(e);
Some(crate::transforms::eval::eval(arena, r))
}
fn expands_to_zero(arena: &mut Arena, e: ExprId) -> bool {
let e = crate::transforms::eval::eval(arena, e);
let e = crate::transforms::expand::expand(arena, e);
let e = crate::transforms::eval::eval(arena, e);
arena.is_zero_structural(e)
}
fn try_linear_cc_symbolic(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
) -> Option<OdeResult> {
let chain = derivative_chain(arena, func, var);
let children: Vec<ExprId> = match arena.node(expr).clone() {
ExprNode::Add(c) => c.to_vec(),
_ => vec![expr],
};
let mut coeffs: Vec<Vec<ExprId>> = vec![Vec::new(); MAX_ODE_ORDER + 1];
let mut forcing = Vec::new();
for child in children {
if let Some((k, c)) = symbolic_linear_term(arena, child, &chain, var, func_sym) {
coeffs[k].push(c);
} else if !contains_sym(arena, child, func_sym) {
forcing.push(child);
} else {
return None;
}
}
let mut a: Vec<ExprId> = Vec::with_capacity(coeffs.len());
for cs in &coeffs {
let s = match cs.len() {
0 => arena.zero,
1 => cs[0],
_ => arena.add(cs),
};
a.push(crate::transforms::eval::eval(arena, s));
}
while a.len() > 1 && a.last().is_some_and(|&c| arena.is_zero_structural(c)) {
a.pop();
}
let n = a.len() - 1;
if n < 2 || a.iter().all(|&c| arena.as_num(c).is_some()) {
return None;
}
if !forcing.is_empty() && n != 2 {
return None;
}
let r = arena.symbol("__r_sym");
let mut terms = Vec::with_capacity(n + 1);
for (k, &c) in a.iter().enumerate() {
let rk = match k {
0 => arena.one,
1 => r,
_ => {
let ke = arena.int(k as i64);
arena.pow(r, ke)
}
};
terms.push(arena.mul(&[c, rk]));
}
let p = arena.add(&terms);
let p = crate::transforms::eval::eval(arena, p);
let (_, factors) = crate::simplify::factor::factor_list(arena, p, Some(r));
let mut basis = Vec::with_capacity(n);
let mut exps: Vec<(ExprId, usize)> = Vec::new();
for (f, m) in factors {
if !crate::base::walk::contains(arena, f, r) {
continue; }
let coeffs = crate::poly::polybridge::poly_coefficients(arena, f, r);
let roots: Vec<ExprId> = match coeffs.as_deref() {
Some([c0, c1]) => {
let q = arena.div(*c0, *c1);
vec![arena.neg(q)]
}
_ => crate::transforms::solve::solve(arena, f, r)
.into_iter()
.map(|s| s.value)
.collect(),
};
if roots.is_empty() {
return None;
}
for r0 in roots {
let r0 = crate::transforms::eval::eval(arena, r0);
if crate::base::walk::contains(arena, r0, var) || exps.iter().any(|&(e, _)| e == r0) {
return None;
}
let mut d = f;
let mut mult = 0usize;
while mult < n {
let at = crate::transforms::subs::subs(arena, d, r, r0);
if !expands_to_zero(arena, at) {
break;
}
mult += 1;
d = crate::transforms::diff::diff(arena, d, r);
}
if mult == 0 {
return None;
}
let mult = mult * m as usize;
exps.push((r0, mult));
for j in 0..mult {
basis.push(mode_real(arena, r0, j, var));
}
}
}
if basis.len() != n {
return None;
}
let homogeneous = real_symbolic_basis(arena, &exps, var).unwrap_or_else(|| basis.clone());
let mut constants = Vec::with_capacity(n);
let mut sum = Vec::with_capacity(n + 1);
for (i, &b) in homogeneous.iter().enumerate() {
let c = arena.symbol(&format!("C{}", i + 1));
constants.push(c);
sum.push(arena.mul(&[c, b]));
}
if !forcing.is_empty() {
let (y1, y2) = (basis[0], basis[1]);
let f = arena.add(&forcing);
let nf = arena.neg(f);
let g = arena.div(nf, a[2]);
let g = crate::transforms::eval::eval(arena, g);
let w = if exps.len() == 2 {
let d = arena.sub(exps[1].0, exps[0].0);
let s = arena.add(&[exps[0].0, exps[1].0]);
let sx = arena.mul(&[s, var]);
let e = arena.exp(sx);
arena.mul(&[d, e])
} else {
let two = arena.int(2);
let sx = arena.mul(&[two, exps[0].0, var]);
arena.exp(sx)
};
let w = crate::transforms::eval::eval(arena, w);
let i1 = {
let num = arena.mul(&[y2, g]);
let q = arena.div(num, w);
let q = crate::transforms::eval::eval(arena, q);
integrate_forms(arena, q, var)
};
if integral_failed(arena, i1) {
return None;
}
let i2 = {
let num = arena.mul(&[y1, g]);
let q = arena.div(num, w);
let q = crate::transforms::eval::eval(arena, q);
integrate_forms(arena, q, var)
};
if integral_failed(arena, i2) {
return None;
}
let t1 = arena.mul(&[y1, i1]);
let t1 = arena.neg(t1);
let t2 = arena.mul(&[y2, i2]);
let yp = arena.add(&[t1, t2]);
let yp = generic_branches(arena, yp, var)?;
let yp = expand_combine(arena, yp);
let yp = crate::poly::polybridge::together(arena, yp);
let yp = crate::transforms::eval::eval(arena, yp);
if crate::base::walk::has_unevaluated(arena, yp) {
return None;
}
sum.push(yp);
}
let solution = arena.add(&sum);
let solution = crate::transforms::eval::eval(arena, solution);
if !forcing.is_empty() && !checkodesol(arena, expr, solution, func, var) {
return None;
}
Some(OdeResult {
solution,
constants,
branches: Vec::new(),
})
}
fn real_symbolic_basis(
arena: &mut Arena,
exps: &[(ExprId, usize)],
var: ExprId,
) -> Option<Vec<ExprId>> {
let mut reals = crate::base::assumptions::AssumptionCache::new();
let parts: Vec<(ExprId, ExprId)> = exps
.iter()
.map(|&(r, _)| {
let post = crate::base::walk::post_order_ids(arena, r);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
let half = Q::new(1.into(), 2.into());
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match *arena.node(rebuilt) {
ExprNode::Pow(b, e) if arena.as_num(e) == Some(&half) => {
let nb = arena.neg(b);
let nb = crate::transforms::eval::eval(arena, nb);
if reals.query(arena, nb, crate::base::assumptions::Props::POSITIVE)
== Some(true)
{
let s = arena.pow(nb, e);
let i_unit = arena.i_unit;
arena.mul(&[i_unit, s])
} else {
rebuilt
}
}
_ => rebuilt,
};
cache.insert(id, new);
}
let r = cache.get(&r).copied().unwrap_or(r);
let r = crate::transforms::eval::eval(arena, r);
let (re, im) = arena.as_real_imag_expr(r);
(
crate::transforms::eval::eval(arena, re),
crate::transforms::eval::eval(arena, im),
)
})
.collect();
let mut used = vec![false; exps.len()];
let mut out = Vec::new();
let mut changed = false;
for i in 0..exps.len() {
if used[i] {
continue;
}
used[i] = true;
let (r, mult) = exps[i];
let (re, im) = parts[i];
let decided = !arena.is_zero_structural(im)
&& real_radical_form(arena, re)
&& real_radical_form(arena, im);
let partner = decided
.then(|| {
(0..exps.len()).find(|&k| {
if used[k] || exps[k].1 != mult {
return false;
}
let (re2, im2) = parts[k];
let dre = arena.sub(re, re2);
let sim = arena.add(&[im, im2]);
expands_to_zero(arena, dre) && expands_to_zero(arena, sim)
})
})
.flatten();
match partner {
Some(k) => {
used[k] = true;
changed = true;
let beta = if arena.as_num(im).is_some_and(|q| q.is_negative())
|| matches!(arena.node(im), ExprNode::Neg(_))
{
parts[k].1
} else {
im
};
for j in 0..mult {
let (c, s) = mode_complex(arena, re, beta, j, var);
out.push(c);
out.push(s);
}
}
None => {
for j in 0..mult {
out.push(mode_real(arena, r, j, var));
}
}
}
}
changed.then_some(out)
}
fn clairaut_f(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
p: ExprId,
) -> Option<ExprId> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if !expr_contains(arena, expr, dy_dx) || expr_contains(arena, expr, d2y_dx2) {
return None;
}
let g = crate::transforms::subs::subs(arena, expr, dy_dx, p);
let coeffs = crate::transforms::solve::symbolic_poly_coeffs(arena, g, func)?;
if coeffs.len() != 2 {
return None;
}
let c_y = coeffs[1];
let c_0 = coeffs[0];
let p_sym = match arena.node(p) {
ExprNode::Symbol(s) => *s,
_ => return None,
};
if contains_sym(arena, c_y, var_sym) || contains_sym(arena, c_y, p_sym) {
return None;
}
let ratio = arena.div(c_0, c_y);
let neg_ratio = arena.neg(ratio);
let xp = arena.mul(&[var, p]);
let f = arena.sub(neg_ratio, xp);
let f = crate::transforms::eval::eval(arena, f);
let f = crate::transforms::expand::expand(arena, f);
let f = crate::transforms::eval::eval(arena, f);
if contains_sym(arena, f, var_sym) || contains_sym(arena, f, func_sym) {
return None;
}
Some(f)
}
fn try_clairaut(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
var_sym: SymbolId,
) -> Option<OdeResult> {
let p = arena.symbol("__clairaut_p");
let f = clairaut_f(arena, expr, func, var, func_sym, var_sym, p)?;
tracing::debug!("ode: Clairaut form recognised");
let c1 = arena.symbol("C1");
let f_c = crate::transforms::subs::subs(arena, f, p, c1);
let cx = arena.mul(&[c1, var]);
let solution = arena.add(&[cx, f_c]);
let solution = crate::transforms::eval::eval(arena, solution);
Some(OdeResult {
solution,
constants: vec![c1],
branches: Vec::new(),
})
}
fn riccati_coeffs(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
) -> Option<(ExprId, ExprId, ExprId)> {
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
if expr_contains(arena, expr, d2y_dx2) {
return None;
}
let (rest, c) = extract_m_n(arena, expr, func, var)?;
if crate::base::walk::contains(arena, c, func) || crate::base::walk::contains(arena, c, var) {
return None;
}
if crate::base::walk::contains(arena, rest, dy_dx) {
return None;
}
let neg_rest = arena.neg(rest);
let rhs = arena.div(neg_rest, c);
let rhs = crate::transforms::eval::eval(arena, rhs);
let coeffs = crate::transforms::solve::symbolic_poly_coeffs(arena, rhs, func)?;
if coeffs.len() != 3 {
return None;
}
let (q0, q1, q2) = (coeffs[0], coeffs[1], coeffs[2]);
if arena.is_zero_structural(q0) || arena.is_zero_structural(q2) {
return None;
}
Some((q0, q1, q2))
}
pub fn solve_riccati(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
particular: ExprId,
) -> Option<OdeResult> {
let (q0, q1, q2) = riccati_coeffs(arena, expr, func, var)?;
if !checkodesol(arena, expr, particular, func, var) {
let yp_prime = crate::transforms::diff::diff(arena, particular, var);
let yp2 = arena.mul(&[particular, particular]);
let t1 = arena.mul(&[q1, particular]);
let t2 = arena.mul(&[q2, yp2]);
let rhs = arena.add(&[q0, t1, t2]);
let residual = arena.sub(yp_prime, rhs);
let residual = crate::transforms::eval::eval(arena, residual);
let residual = crate::transforms::expand::expand(arena, residual);
let residual = crate::transforms::eval::eval(arena, residual);
if !arena.is_zero_structural(residual) {
return None;
}
}
let _ = q0;
let v = arena.symbol("__riccati_v");
let dv = arena.intern(ExprNode::Derivative(v, var));
let two = arena.int(2);
let two_q2_yp = arena.mul(&[two, q2, particular]);
let coeff = arena.add(&[q1, two_q2_yp]);
let coeff = crate::transforms::eval::eval(arena, coeff);
let coeff_v = arena.mul(&[coeff, v]);
let lin = arena.add(&[dv, coeff_v, q2]);
let lin = crate::transforms::eval::eval(arena, lin);
let v_res = dsolve(arena, lin, v, var)?;
if crate::base::walk::contains(arena, v_res.solution, v) {
return None; }
let neg_one = arena.neg_one;
let inv_v = arena.pow(v_res.solution, neg_one);
let solution = arena.add(&[particular, inv_v]);
let solution = crate::transforms::eval::eval(arena, solution);
Some(OdeResult {
solution,
constants: v_res.constants,
branches: Vec::new(),
})
}
#[derive(Clone, Debug, PartialEq)]
struct RatFn {
num: crate::poly::Poly,
den: crate::poly::Poly,
}
impl RatFn {
fn new(num: crate::poly::Poly, den: crate::poly::Poly) -> Option<Self> {
use crate::poly::Poly;
if den.is_zero() {
return None;
}
if num.is_zero() {
return Some(Self {
num: Poly::zero(),
den: Poly::one(),
});
}
let g = Poly::gcd(&num, &den);
let (num, den) = if g.degree().unwrap_or(0) > 0 {
(num.div(&g), den.div(&g))
} else {
(num, den)
};
let inv = Q::one() / den.leading_coeff()?.clone();
Some(Self {
num: num.scale(&inv),
den: den.scale(&inv),
})
}
fn poly(p: crate::poly::Poly) -> Self {
Self {
num: p,
den: crate::poly::Poly::one(),
}
}
fn constant(c: Q) -> Self {
Self::poly(crate::poly::Poly::constant(c))
}
fn is_zero(&self) -> bool {
self.num.is_zero()
}
fn add(&self, o: &Self) -> Option<Self> {
Self::new(
self.num.mul(&o.den).add(&o.num.mul(&self.den)),
self.den.mul(&o.den),
)
}
fn sub(&self, o: &Self) -> Option<Self> {
self.add(&o.scale(&-Q::one()))
}
fn mul(&self, o: &Self) -> Option<Self> {
Self::new(self.num.mul(&o.num), self.den.mul(&o.den))
}
fn div(&self, o: &Self) -> Option<Self> {
if o.is_zero() {
return None;
}
Self::new(self.num.mul(&o.den), self.den.mul(&o.num))
}
fn scale(&self, c: &Q) -> Self {
Self::new(self.num.scale(c), self.den.clone()).unwrap_or_else(|| self.clone())
}
fn deriv(&self) -> Option<Self> {
let n = self
.num
.derivative()
.mul(&self.den)
.sub(&self.num.mul(&self.den.derivative()));
Self::new(n, self.den.mul(&self.den))
}
fn val_inf(&self) -> i64 {
let d = |p: &crate::poly::Poly| i64::try_from(p.degree().unwrap_or(0)).unwrap_or(0);
if self.is_zero() {
return i64::MAX / 4;
}
d(&self.den) - d(&self.num)
}
fn to_expr(&self, arena: &mut Arena, var: ExprId) -> ExprId {
let n = crate::poly::polybridge::poly_to_expr(arena, &self.num, var);
if self.den.is_constant() {
return n;
}
let d = crate::poly::polybridge::poly_to_expr(arena, &self.den, var);
arena.div(n, d)
}
fn from_expr(arena: &mut Arena, e: ExprId, var: ExprId) -> Option<Self> {
let t = crate::poly::polybridge::together(arena, e);
let (n, d) = crate::poly::polybridge::as_numer_denom(arena, t);
let n = crate::poly::polybridge::expr_to_poly(arena, n, var)?;
let d = crate::poly::polybridge::expr_to_poly(arena, d, var)?;
Self::new(n, d)
}
}
fn power_series(p: &crate::poly::Poly, q: &crate::poly::Poly, count: usize) -> Option<Vec<Q>> {
use num_traits::Zero;
let q0 = q.coeff(0);
if q0.is_zero() {
return None;
}
let qd = q.degree().unwrap_or(0);
let mut g: Vec<Q> = Vec::with_capacity(count);
for k in 0..count {
let mut s = p.coeff(k);
for j in 1..=k.min(qd) {
s -= q.coeff(j) * &g[k - j];
}
g.push(s / &q0);
}
Some(g)
}
fn laurent_at_pole(f: &RatFn, p: &Q, m: usize, count: usize) -> Option<Vec<Q>> {
use num_traits::Zero;
let num = f.num.taylor_shift(p);
let den = f.den.taylor_shift(p);
if (0..m).any(|k| !den.coeff(k).is_zero()) {
return None;
}
let den1 = crate::poly::Poly::from_coeffs(den.coeffs().get(m..)?.to_vec());
power_series(&num, &den1, count)
}
fn laurent_at_infinity(f: &RatFn, count: usize) -> Option<Vec<Q>> {
let rev = |p: &crate::poly::Poly| {
let mut c = p.coeffs().to_vec();
c.reverse();
crate::poly::Poly::from_coeffs(c)
};
power_series(&rev(&f.num), &rev(&f.den), count)
}
fn rational_sqrt(q: &Q) -> Option<Q> {
if q.is_negative() {
return None;
}
let (n, d) = (q.numer(), q.denom());
let (rn, rd) = (n.sqrt(), d.sqrt());
(&rn * &rn == *n && &rd * &rd == *d).then(|| Q::new(rn, rd))
}
const MAX_RICCATI_POLES: usize = 6;
const MAX_RICCATI_AUX_DEGREE: i64 = 24;
fn rational_riccati_particulars(a: &RatFn) -> Option<Vec<RatFn>> {
use crate::poly::Poly;
use num_traits::{ToPrimitive, Zero};
let one = Q::one();
let two = Q::from_integer(2.into());
let mut poles: Vec<(Q, usize)> = Vec::new();
if !a.den.is_constant() {
let (_, factors) = a.den.factor_over_z();
for (f, mult) in factors {
if f.degree() != Some(1) {
return None;
}
poles.push((-(f.coeff(0) / f.coeff(1)), usize::try_from(mult).ok()?));
}
}
if poles.len() > MAX_RICCATI_POLES {
return None;
}
let val_inf = a.val_inf();
if poles.iter().any(|&(_, m)| m != 1 && !m.is_multiple_of(2))
|| !(val_inf >= 2 || (val_inf <= 0 && val_inf % 2 == 0))
{
return Some(Vec::new());
}
let mut c_choices: Vec<Vec<Vec<Q>>> = Vec::with_capacity(poles.len());
for (p, m) in &poles {
let mut choices = Vec::new();
if *m == 1 {
choices.push(vec![one.clone()]);
} else if *m == 2 {
let ser = laurent_at_pole(a, p, 2, 1)?;
let disc = &one + Q::from_integer(4.into()) * &ser[0];
if disc.is_zero() {
choices.push(vec![&one / &two]);
} else if let Some(s) = rational_sqrt(&disc) {
choices.push(vec![(&one + &s) / &two]);
choices.push(vec![(&one - &s) / &two]);
}
} else {
let ri = m / 2;
let raw = laurent_at_pole(a, p, *m, *m + 1)?;
let ser = |key: usize| raw.get(*m - key).cloned().unwrap_or_else(Q::zero);
let Some(lead) = rational_sqrt(&ser(2 * ri)) else {
c_choices.push(Vec::new());
continue;
};
for sign in [one.clone(), -one.clone()] {
let mut c = vec![Q::zero(); ri];
c[ri - 1] = &sign * &lead;
if c[ri - 1].is_zero() {
continue;
}
let mut sm = Q::zero();
let mut s_last = ri;
for s in (1..ri).rev() {
sm = Q::zero();
for j in s + 1..ri {
sm += &c[j - 1] * &c[ri + s - j - 1];
}
if s != 1 {
c[s - 1] = (ser(ri + s) - &sm) / (&two * &c[ri - 1]);
}
s_last = s;
}
let ri_q = Q::from_integer(ri.into());
c[0] = (ser(ri + s_last) - &sm - &ri_q * &c[ri - 1]) / (&two * &c[ri - 1]);
if !choices.contains(&c) {
choices.push(c);
}
}
}
c_choices.push(choices);
}
let mut d_choices: Vec<(Vec<Q>, Q)> = Vec::new();
if val_inf < 0 {
let n = usize::try_from(-val_inf / 2).ok()?;
let r = laurent_at_infinity(a, n + 2)?;
let ser = |j: i64| -> Q {
let i = 2 * i64::try_from(n).unwrap_or(0) - j;
usize::try_from(i)
.ok()
.and_then(|i| r.get(i).cloned())
.unwrap_or_else(Q::zero)
};
let n_i = i64::try_from(n).ok()?;
if let Some(lead) = rational_sqrt(&ser(2 * n_i)) {
for sign in [one.clone(), -one.clone()] {
let mut d = vec![Q::zero(); n + 1];
d[n] = &sign * &lead;
if d[n].is_zero() {
continue;
}
let mut sm = Q::zero();
for s in (-1..n_i).rev() {
sm = Q::zero();
for j in (s + 1)..n_i {
let (ju, ku) =
(usize::try_from(j).ok()?, usize::try_from(n_i + s - j).ok()?);
sm += &d[ju] * &d[ku];
}
if s >= 0 {
let su = usize::try_from(s).ok()?;
d[su] = (ser(n_i + s) - &sm) / (&two * &d[n]);
}
}
let n_q = Q::from_integer(n.into());
let dm1 = (ser(n_i - 1) - &n_q * &d[n] - &sm) / (&two * &d[n]);
if !d_choices.iter().any(|(v, w)| *v == d && *w == dm1) {
d_choices.push((d, dm1));
}
}
}
} else if val_inf == 0 {
let r = laurent_at_infinity(a, 2)?;
if let Some(d0) = rational_sqrt(&r[0])
&& !d0.is_zero()
{
for sign in [one.clone(), -one.clone()] {
let d0s = &sign * &d0;
let dm1 = &r[1] / (&two * &d0s);
d_choices.push((vec![d0s], dm1));
}
}
} else {
let s_inf = if val_inf == 2 {
laurent_at_infinity(a, 1)?[0].clone()
} else {
Q::zero()
};
let disc = &one + Q::from_integer(4.into()) * &s_inf;
if disc.is_zero() {
d_choices.push((Vec::new(), &one / &two));
} else if let Some(s) = rational_sqrt(&disc) {
d_choices.push((Vec::new(), (&one + &s) / &two));
d_choices.push((Vec::new(), (&one - &s) / &two));
}
}
let mut found: Vec<RatFn> = Vec::new();
let mut idx = vec![0usize; c_choices.len()];
if c_choices.iter().any(Vec::is_empty) || d_choices.is_empty() {
return Some(found);
}
loop {
for (d, dm1) in &d_choices {
let mut m = dm1.clone();
for (k, ch) in c_choices.iter().enumerate() {
m -= &ch[idx[k]][0];
}
if m.is_integer() && !m.is_negative() {
let m = m.to_integer().to_i64()?;
if m <= MAX_RICCATI_AUX_DEGREE {
let mut ybar = RatFn::poly(Poly::from_coeffs(d.clone()));
for (k, (p, _)) in poles.iter().enumerate() {
let lin = Poly::from_coeffs(vec![-p.clone(), one.clone()]);
let mut pw = Poly::one();
for cij in &c_choices[k][idx[k]] {
pw = pw.mul(&lin);
let term = RatFn::new(Poly::constant(cij.clone()), pw.clone())?;
ybar = ybar.add(&term)?;
}
}
if let Some(y) = riccati_aux_solution(a, &ybar, usize::try_from(m).ok()?)
&& !found.contains(&y)
{
found.push(y);
}
}
}
}
let mut k = 0;
loop {
if k == idx.len() {
return Some(found);
}
idx[k] += 1;
if idx[k] < c_choices[k].len() {
break;
}
idx[k] = 0;
k += 1;
}
}
}
fn riccati_aux_solution(a: &RatFn, ybar: &RatFn, m: usize) -> Option<RatFn> {
use crate::poly::Poly;
use num_traits::Zero;
let a0 = ybar.deriv()?.add(&ybar.mul(ybar)?)?.sub(a)?;
let a1 = ybar.scale(&Q::from_integer(2.into()));
let lcm = {
let g = Poly::gcd(&a0.den, &a1.den);
a0.den.mul(&a1.den).div(&g)
};
let b0 = a0.num.mul(&lcm.div(&a0.den));
let b1 = a1.num.mul(&lcm.div(&a1.den));
let apply = |p: &Poly| -> Poly {
let d1 = p.derivative();
let d2 = d1.derivative();
lcm.mul(&d2).add(&b1.mul(&d1)).add(&b0.mul(p))
};
let mono = |k: usize| Poly::monomial(Q::one(), k);
let rhs = apply(&mono(m)).neg();
let cols: Vec<Poly> = (0..m).map(|k| apply(&mono(k))).collect();
let rows = cols
.iter()
.map(|c| c.degree().unwrap_or(0))
.chain(std::iter::once(rhs.degree().unwrap_or(0)))
.max()
.unwrap_or(0)
+ 1;
let mut mat: Vec<Vec<Q>> = (0..rows)
.map(|r| {
let mut row: Vec<Q> = cols.iter().map(|c| c.coeff(r)).collect();
row.push(rhs.coeff(r));
row
})
.collect();
let mut pivots: Vec<(usize, usize)> = Vec::new();
let mut row = 0;
for col in 0..m {
let Some(pr) = (row..rows).find(|&r| !mat[r][col].is_zero()) else {
continue;
};
mat.swap(row, pr);
let inv = Q::one() / mat[row][col].clone();
for v in mat[row].iter_mut() {
*v *= &inv;
}
let pivot_row = mat[row].clone();
for (r, target) in mat.iter_mut().enumerate() {
if r != row && !target[col].is_zero() {
let f = target[col].clone();
for (t, p) in target.iter_mut().zip(&pivot_row) {
*t -= &f * p;
}
}
}
pivots.push((row, col));
row += 1;
}
if (row..rows).any(|r| !mat[r][m].is_zero()) {
return None;
}
let mut coeffs = vec![Q::zero(); m + 1];
coeffs[m] = Q::one();
for (r, c) in pivots {
coeffs[c] = mat[r][m].clone();
}
let p = Poly::from_coeffs(coeffs);
let y = ybar.add(&RatFn::new(p.derivative(), p)?)?;
y.deriv()?.add(&y.mul(&y)?)?.sub(a)?.is_zero().then_some(y)
}
fn try_rational_riccati(
arena: &mut Arena,
expr: ExprId,
func: ExprId,
var: ExprId,
func_sym: SymbolId,
_var_sym: SymbolId,
) -> Option<OdeResult> {
let rhs = first_order_rhs(arena, expr, func, var)?;
let coeffs = crate::transforms::solve::symbolic_poly_coeffs(arena, rhs, func)?;
if coeffs.len() != 3 || coeffs.iter().any(|&c| contains_sym(arena, c, func_sym)) {
return None;
}
for &c in &coeffs {
if crate::base::walk::free_symbols(arena, c)
.iter()
.any(|&s| s != var)
{
return None;
}
}
let b0 = RatFn::from_expr(arena, coeffs[0], var)?;
let b1 = RatFn::from_expr(arena, coeffs[1], var)?;
let b2 = RatFn::from_expr(arena, coeffs[2], var)?;
if b0.is_zero() || b2.is_zero() {
return None;
}
let q = |n: i64, d: i64| Q::new(n.into(), d.into());
let b1p = b1.deriv()?;
let b2p = b2.deriv()?;
let b2pp = b2p.deriv()?;
let a = b0
.mul(&b2)?
.scale(&q(-1, 1))
.add(&b1.mul(&b1)?.scale(&q(1, 4)))?
.sub(&b1p.scale(&q(1, 2)))?
.add(&b2p.mul(&b2p)?.div(&b2.mul(&b2)?)?.scale(&q(3, 4)))?
.add(&b1.mul(&b2p)?.div(&b2)?.scale(&q(1, 2)))?
.sub(&b2pp.div(&b2)?.scale(&q(1, 2)))?;
let c1 = arena.symbol("C1");
let v = if a.is_zero() {
let s = arena.add(&[var, c1]);
let m1 = arena.int(-1);
arena.pow(s, m1)
} else {
let parts = rational_riccati_particulars(&a)?;
let parts: Vec<ExprId> = parts.iter().map(|p| p.to_expr(arena, var)).collect();
match parts.as_slice() {
[] => return None,
[v1] => {
let two = arena.int(2);
let tv = arena.mul(&[two, *v1]);
let f = integrate_forms(arena, tv, var);
if integral_failed(arena, f) {
return None;
}
let ef = exp_of_log_sum(arena, f);
let m1 = arena.int(-1);
let emf = arena.pow(ef, m1);
let emf = crate::transforms::eval::eval(arena, emf);
let g = integrate_forms(arena, emf, var);
if integral_failed(arena, g) {
return None;
}
let inner = arena.add(&[c1, g]);
let z = arena.mul(&[ef, inner]);
let inv_z = arena.pow(z, m1);
arena.add(&[*v1, inv_z])
}
[v1, v2, ..] => {
let neg_v1 = arena.neg(*v1);
let diff = arena.add(&[*v2, neg_v1]);
let diff = crate::transforms::eval::eval(arena, diff);
let f = integrate_forms(arena, diff, var);
if integral_failed(arena, f) {
return None;
}
let ef = exp_of_log_sum(arena, f);
let u = arena.mul(&[c1, ef]);
let v2u = arena.mul(&[*v2, u]);
let num = arena.add(&[v2u, neg_v1]);
let m1 = arena.int(-1);
let den = arena.add(&[u, m1]);
arena.div(num, den)
}
}
};
let bp = b2p
.div(&b2.mul(&b2)?)?
.scale(&q(-1, 2))
.sub(&b1.div(&b2)?.scale(&q(1, 2)))?;
let inv_b2 = RatFn::constant(Q::one()).div(&b2)?.scale(&q(-1, 1));
let inv_b2 = inv_b2.to_expr(arena, var);
let bp = bp.to_expr(arena, var);
let vy = arena.mul(&[inv_b2, v]);
let sol = arena.add(&[vy, bp]);
let sol = crate::transforms::eval::eval(arena, sol);
let sol = crate::simplify::ratsimp::ratsimp(arena, sol);
let sol = crate::transforms::eval::eval(arena, sol);
if !crate::base::walk::contains(arena, sol, c1)
|| crate::base::walk::has_unevaluated(arena, sol)
|| !family_verifies(arena, expr, sol, func, var, c1)
{
return None;
}
Some(OdeResult {
solution: sol,
constants: vec![c1],
branches: Vec::new(),
})
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum OdeType {
SimpleSeparable,
FullSeparable,
FirstOrderLinearCC,
FirstOrderLinearVC,
ExactFirstOrder,
SecondOrderLinearCCHomogeneous,
SecondOrderLinearCCNonHomogeneous,
Bernoulli,
EulerCauchy,
VariationOfParameters,
HomogeneousCoefficient,
NthOrderReducible,
NthOrderLinearConstCoeff,
IntegratingFactor,
Clairaut,
Riccati,
Unknown,
}
pub fn classify_ode(arena: &mut Arena, expr: ExprId, func: ExprId, var: ExprId) -> OdeType {
let var_sym = match arena.node(var) {
ExprNode::Symbol(sid) => *sid,
_ => return OdeType::Unknown,
};
let func_sym = match arena.node(func) {
ExprNode::Symbol(sid) => *sid,
_ => return OdeType::Unknown,
};
let dy_dx = arena.intern(ExprNode::Derivative(func, var));
let d2y_dx2 = arena.intern(ExprNode::Derivative(dy_dx, var));
let d3y_dx3 = arena.intern(ExprNode::Derivative(d2y_dx2, var));
if expr_contains(arena, expr, d3y_dx3) {
if let Some(ode) = extract_linear_cc(arena, expr, func, var, func_sym)
&& ode.coeffs.len() >= 4
{
return OdeType::NthOrderLinearConstCoeff;
}
return OdeType::Unknown;
}
if expr_contains(arena, expr, d2y_dx2) {
if let ExprNode::Add(ref children) = arena.node(expr).clone() {
let mut all_const_coeff = true;
let mut has_y2 = false;
let mut forcing: Vec<ExprId> = Vec::new();
for &child in children {
let (_coeff, term) = arena.as_coeff_term(child);
if term == d2y_dx2 {
has_y2 = true;
} else if term == dy_dx {
} else if term == func {
} else if !contains_sym(arena, child, func_sym) {
forcing.push(child);
} else {
all_const_coeff = false;
break;
}
}
if has_y2 && all_const_coeff {
if forcing.is_empty() {
return OdeType::SecondOrderLinearCCHomogeneous;
}
let all_uc = forcing
.iter()
.all(|&t| parse_forcing_term(arena, t, var).is_some());
if all_uc {
return OdeType::SecondOrderLinearCCNonHomogeneous;
}
if try_variation_of_parameters(arena, expr, func, var, func_sym, var_sym).is_some()
{
return OdeType::VariationOfParameters;
}
return OdeType::Unknown;
}
}
if try_euler_cauchy(arena, expr, func, var, func_sym, var_sym).is_some() {
return OdeType::EulerCauchy;
}
{
let d2_ph = arena.symbol("__d2_cls");
let d1_ph = arena.symbol("__d1_cls");
let stripped = crate::transforms::subs::subs(arena, expr, d2y_dx2, d2_ph);
let stripped = crate::transforms::subs::subs(arena, stripped, dy_dx, d1_ph);
if !contains_sym(arena, stripped, var_sym) {
return OdeType::NthOrderReducible;
}
}
return OdeType::Unknown;
}
if expr_contains(arena, expr, dy_dx) {
if !contains_sym_outside_deriv(arena, expr, func_sym, dy_dx) {
return OdeType::SimpleSeparable;
}
if let ExprNode::Add(ref children) = arena.node(expr).clone() {
let mut is_linear = true;
let mut has_var_coeff = false;
for &child in children {
let (coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
} else if term == func {
} else if contains_sym(arena, child, func_sym) {
if let Some(px) = extract_coeff_of_func(arena, child, func, func_sym, var_sym) {
let _ = coeff; if contains_sym(arena, px, var_sym) {
has_var_coeff = true;
}
} else {
is_linear = false;
break;
}
}
}
if is_linear {
if has_var_coeff {
return OdeType::FirstOrderLinearVC;
}
return OdeType::FirstOrderLinearCC;
}
}
if let Some((m_ex, n_ex)) = extract_m_n(arena, expr, func, var)
&& (contains_sym(arena, m_ex, func_sym) || contains_sym(arena, n_ex, func_sym))
{
let dm_dy = crate::transforms::diff::diff(arena, m_ex, func);
let dn_dx = crate::transforms::diff::diff(arena, n_ex, var);
let check = arena.sub(dm_dy, dn_dx);
let check = crate::transforms::eval::eval(arena, check);
let check = crate::transforms::expand::expand(arena, check);
let check = crate::transforms::eval::eval(arena, check);
if check == arena.zero {
return OdeType::ExactFirstOrder;
}
}
{
let p = arena.symbol("__clairaut_p_cls");
if clairaut_f(arena, expr, func, var, func_sym, var_sym, p).is_some() {
return OdeType::Clairaut;
}
}
if let ExprNode::Add(ref bn_children) = arena.node(expr).clone() {
let mut bn_has_dy = false;
let mut bn_has_yn = false;
let mut bn_ok = true;
let mut bn_has_free = false;
for &child in bn_children {
let (_, term) = arena.as_coeff_term(child);
if term == dy_dx {
bn_has_dy = true;
} else if !contains_sym(arena, child, func_sym) {
bn_has_free = true;
} else if extract_coeff_of_func(arena, child, func, func_sym, var_sym).is_some() {
} else if extract_bernoulli_term(arena, child, func, func_sym, var_sym).is_some() {
bn_has_yn = true;
} else {
bn_ok = false;
break;
}
}
if bn_has_dy && bn_has_yn && bn_ok && !bn_has_free {
return OdeType::Bernoulli;
}
}
if let ExprNode::Add(ref hc_children) = arena.node(expr).clone() {
let mut hc_has_dy = false;
let mut hc_other: Vec<ExprId> = Vec::new();
for &child in hc_children {
let (_, term) = arena.as_coeff_term(child);
if term == dy_dx {
hc_has_dy = true;
} else {
hc_other.push(child);
}
}
if hc_has_dy && !hc_other.is_empty() {
let hc_rhs = if hc_other.len() == 1 {
arena.neg(hc_other[0])
} else {
let s = arena.add(&hc_other);
arena.neg(s)
};
if contains_sym(arena, hc_rhs, func_sym) && contains_sym(arena, hc_rhs, var_sym) {
let v_cls = arena.symbol("__v_cls");
let vx = arena.mul(&[v_cls, var]);
let sub = crate::transforms::subs::subs(arena, hc_rhs, func, vx);
let sub = crate::transforms::eval::eval(arena, sub);
let sub = crate::transforms::expand::expand(arena, sub);
let sub = crate::transforms::eval::eval(arena, sub);
let sub = arena.cancel_expr(sub, var);
let sub = crate::transforms::eval::eval(arena, sub);
if !contains_sym(arena, sub, var_sym) {
return OdeType::HomogeneousCoefficient;
}
}
}
}
if let ExprNode::Add(ref children) = arena.node(expr).clone() {
let mut has_deriv = false;
let mut other_terms: Vec<ExprId> = Vec::new();
for &child in children {
let (_coeff, term) = arena.as_coeff_term(child);
if term == dy_dx {
has_deriv = true;
} else {
other_terms.push(child);
}
}
if has_deriv && !other_terms.is_empty() {
let rhs = if other_terms.len() == 1 {
arena.neg(other_terms[0])
} else {
let sum = arena.add(&other_terms);
arena.neg(sum)
};
let factors = collect_mul_factors(arena, rhs);
let mut can_separate = true;
for &factor in &factors {
let has_x = contains_sym(arena, factor, var_sym);
let has_y = contains_sym(arena, factor, func_sym);
if has_x && has_y {
can_separate = false;
break;
}
}
if can_separate {
return OdeType::FullSeparable;
}
}
}
if try_integrating_factor_ode(arena, expr, func, var, func_sym, var_sym).is_some() {
return OdeType::IntegratingFactor;
}
if riccati_coeffs(arena, expr, func, var).is_some() {
return OdeType::Riccati;
}
return OdeType::Unknown;
}
OdeType::Unknown
}
fn expr_contains(arena: &Arena, haystack: ExprId, needle: ExprId) -> bool {
crate::base::walk::contains(arena, haystack, needle)
}
fn contains_sym_outside_deriv(
arena: &Arena,
expr: ExprId,
sym: SymbolId,
deriv_node: ExprId,
) -> bool {
let mut stack = vec![expr];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(id) = stack.pop() {
if id == deriv_node || !seen.insert(id) {
continue;
}
let node = arena.node(id);
if let ExprNode::Symbol(s) = node {
if *s == sym {
return true;
}
continue;
}
node.for_each_child(|c| stack.push(c));
}
false
}
pub fn checkodesol(
arena: &mut Arena,
ode_expr: ExprId,
solution: ExprId,
func: ExprId,
var: ExprId,
) -> bool {
let chain = derivative_chain(arena, func, var);
let max_order = (1..chain.len())
.rev()
.find(|&k| expr_contains(arena, ode_expr, chain[k]))
.unwrap_or(0);
let mut sol_derivs = vec![solution];
for k in 1..=max_order {
let d = crate::transforms::diff::diff(arena, sol_derivs[k - 1], var);
sol_derivs.push(d);
}
let mut result = ode_expr;
for k in (0..=max_order).rev() {
result = crate::transforms::subs::subs(arena, result, chain[k], sol_derivs[k]);
}
result = crate::transforms::eval::eval(arena, result);
result = crate::transforms::expand::expand(arena, result);
result = crate::transforms::eval::eval(arena, result);
if result == arena.zero {
return true;
}
result = crate::transforms::expand::expand(arena, result);
result = crate::transforms::eval::eval(arena, result);
if result == arena.zero {
return true;
}
let simplified = crate::simplify::simplify_engine::unified_simplify(
arena,
result,
&crate::simplify::simplify_engine::SimplifyOpts::default(),
);
if arena.is_zero_structural(simplified.expr) {
return true;
}
radical_normal_form_is_zero(arena, result) || trig_normal_form_is_zero(arena, result)
}
const MAX_RADICAL_TERMS: usize = 64;
type TermRadicals = (Vec<ExprId>, Vec<(ExprId, Q)>);
fn expand_combine(arena: &mut Arena, e: ExprId) -> ExprId {
let e = crate::transforms::expand::expand(arena, e);
let e = crate::transforms::eval::eval(arena, e);
let e = crate::simplify::powsimp::powsimp(arena, e);
crate::transforms::eval::eval(arena, e)
}
fn radical_normal_form_is_zero(arena: &mut Arena, r: ExprId) -> bool {
let r = expand_combine(arena, r);
let terms: Vec<ExprId> = match arena.node(r) {
ExprNode::Add(cs) => cs.to_vec(),
_ => vec![r],
};
if terms.len() > MAX_RADICAL_TERMS {
return false;
}
let mut split: Vec<TermRadicals> = Vec::with_capacity(terms.len());
let mut bases: Vec<ExprId> = Vec::new();
for &t in &terms {
let factors: Vec<ExprId> = match arena.node(t) {
ExprNode::Mul(fs) => fs.to_vec(),
_ => vec![t],
};
let mut rest = Vec::new();
let mut rads = Vec::new();
for f in factors {
match arena.node(f).clone() {
ExprNode::Pow(b, e)
if arena.as_num(e).is_some_and(|q| !q.is_integer())
&& !matches!(arena.node(b), ExprNode::Num(_)) =>
{
let q = arena.as_num(e).cloned().unwrap_or_default();
let nb = expand_combine(arena, b);
if !bases.contains(&nb) {
bases.push(nb);
}
rads.push((nb, q));
}
_ => rest.push(f),
}
}
split.push((rest, rads));
}
if bases.is_empty() || bases.len() > MAX_RADICAL_TERMS {
return false;
}
let exponent = |rads: &[(ExprId, Q)], b: ExprId| -> Q {
rads.iter()
.filter(|(nb, _)| *nb == b)
.map(|(_, q)| q.clone())
.sum()
};
let mut minima = Vec::with_capacity(bases.len());
for &b in &bases {
let mut m: Option<Q> = None;
for (_, rads) in &split {
let e = exponent(rads, b);
m = Some(match m {
Some(v) if v <= e => v,
_ => e,
});
}
minima.push(m.unwrap_or_default());
}
let mut new_terms = Vec::with_capacity(split.len());
for (rest, rads) in &split {
let mut fs = rest.clone();
for (&b, m) in bases.iter().zip(&minima) {
let k = exponent(rads, b) - m;
if !k.is_integer() {
return false;
}
if num_traits::Zero::is_zero(&k) {
continue;
}
let Some(ki) = num_traits::ToPrimitive::to_i64(&k.to_integer()) else {
return false;
};
if ki > 8 {
return false;
}
let ke = arena.int(ki);
fs.push(arena.pow(b, ke));
}
new_terms.push(arena.mul(&fs));
}
let s = arena.add(&new_terms);
let s = expand_combine(arena, s);
arena.is_zero_structural(s) || trig_normal_form_is_zero(arena, s)
}
fn trig_normal_form_is_zero(arena: &mut Arena, r: ExprId) -> bool {
let post = crate::base::walk::post_order_ids(arena, r);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match arena.node(rebuilt).clone() {
ExprNode::Tan(u) => {
let s = arena.sin(u);
let c = arena.cos(u);
arena.div(s, c)
}
_ => rebuilt,
};
cache.insert(id, new);
}
let t = cache.get(&r).copied().unwrap_or(r);
let t = crate::transforms::eval::eval(arena, t);
let t = crate::poly::polybridge::together(arena, t);
let (num, _) = crate::poly::polybridge::as_numer_denom(arena, t);
let mut num = expand_combine(arena, num);
for _ in 0..6 {
if arena.is_zero_structural(num) {
return true;
}
let post = crate::base::walk::post_order_ids(arena, num);
let mut cache: rustc_hash::FxHashMap<ExprId, ExprId> = rustc_hash::FxHashMap::default();
let mut changed = false;
for &id in &post {
let rebuilt = crate::base::walk::rebuild_with_cache(arena, id, &cache);
let new = match arena.node(rebuilt).clone() {
ExprNode::Pow(b, e) => {
let k = arena.as_num(e).and_then(|q| {
q.is_integer()
.then(|| num_traits::ToPrimitive::to_i64(&q.to_integer()))
.flatten()
});
match (arena.node(b).clone(), k) {
(ExprNode::Sin(u), Some(k)) if (2..=64).contains(&k) => {
changed = true;
let c = arena.cos(u);
let two = arena.int(2);
let c2 = arena.pow(c, two);
let one = arena.one;
let one_minus = arena.sub(one, c2);
let half = arena.int(k / 2);
let p = arena.pow(one_minus, half);
if k % 2 == 1 { arena.mul(&[b, p]) } else { p }
}
_ => rebuilt,
}
}
_ => rebuilt,
};
cache.insert(id, new);
}
if !changed {
return false;
}
let n = cache.get(&num).copied().unwrap_or(num);
num = expand_combine(arena, n);
}
arena.is_zero_structural(num)
}
pub fn solve_ode_system(a_matrix: &Matrix, t_var: &Ex) -> Option<Vec<Ex>> {
let n = a_matrix.nrows();
if !a_matrix.is_square() || n == 0 {
return None;
}
for i in 0..n {
for j in 0..n {
if a_matrix.get(i, j).contains(t_var) {
return None;
}
}
}
if ode_system_is_diagonal(a_matrix, n) {
return Some(solve_ode_system_diagonal(a_matrix, t_var, n));
}
if let Some(sol) = solve_ode_system_eigen(a_matrix, t_var, n) {
return Some(sol);
}
solve_ode_system_expm(a_matrix, t_var, n)
}
pub fn solve_ode_system_ivp(
a_matrix: &Matrix,
t_var: &Ex,
x0: &[Ex],
) -> Result<Vec<Ex>, crate::base::errors::SymplexError> {
use crate::base::errors::SymplexError;
let n = a_matrix.nrows();
if !a_matrix.is_square() || n == 0 {
return Err(SymplexError::InvalidArgument {
operation: "solve_ode_system_ivp",
reason: "coefficient matrix must be square and non-empty".into(),
});
}
if x0.len() != n {
return Err(SymplexError::InvalidArgument {
operation: "solve_ode_system_ivp",
reason: format!("expected {n} initial values, got {}", x0.len()),
});
}
let general =
solve_ode_system(a_matrix, t_var).ok_or_else(|| SymplexError::ComputationFailed {
operation: "solve_ode_system_ivp",
reason: "could not solve the homogeneous system".into(),
})?;
let ctx = t_var.context();
let constants: Vec<Ex> = (1..=n).map(|i| ctx.symbol(&format!("C{i}"))).collect();
let zero = ctx.int(0);
let eqs: Vec<Ex> = general
.iter()
.zip(x0)
.map(|(g, v)| (g.subs(t_var, &zero).eval() - v).eval())
.collect();
let sol = crate::api::expr_solve_ext::linsolve(&eqs, &constants).map_err(|e| {
SymplexError::ComputationFailed {
operation: "solve_ode_system_ivp",
reason: format!("could not fit initial values: {e}"),
}
})?;
let pairs = match sol {
crate::api::expr_solve_ext::LinearSolution::Inconsistent => {
return Err(SymplexError::NoSolution {
operation: "solve_ode_system_ivp",
reason: "initial values are inconsistent with the general solution".into(),
});
}
crate::api::expr_solve_ext::LinearSolution::Unique(p) => p,
crate::api::expr_solve_ext::LinearSolution::Parametric { solution, free } => solution
.into_iter()
.filter(|(c, _)| !free.contains(c))
.collect(),
};
Ok(general
.iter()
.map(|g| {
let mut e = g.clone();
for (c, v) in &pairs {
e = e.subs(c, v);
}
e.eval().simplify()
})
.collect())
}
pub fn solve_ode_system_nonhomogeneous(
a_matrix: &Matrix,
b_vec: &[Ex],
t_var: &Ex,
) -> Option<Vec<Ex>> {
let n = a_matrix.nrows();
if !a_matrix.is_square() || n == 0 || b_vec.len() != n {
return None;
}
let x_h = solve_ode_system(a_matrix, t_var)?;
let neg_a = a_matrix.scale(&t_var.context().int(-1));
let exp_neg_at = neg_a.matrix_exp_t(t_var).ok()?;
let b_col = Matrix::col_vector(b_vec.to_vec()).ok()?;
let integrand_matrix = exp_neg_at.matmul(&b_col).ok()?.eval();
let mut integrated = Vec::with_capacity(n);
for i in 0..n {
integrated.push(integrand_matrix.get(i, 0).integrate(t_var).eval());
}
let integrated_col = Matrix::col_vector(integrated).ok()?;
let exp_at = a_matrix.matrix_exp_t(t_var).ok()?;
let particular = exp_at.matmul(&integrated_col).ok()?.eval();
let mut solution = Vec::with_capacity(n);
for (i, x_h_i) in x_h.iter().enumerate() {
let xi: Ex = x_h_i + particular.get(i, 0);
solution.push(xi.eval());
}
Some(solution)
}
pub fn classify_ode_system_is_constant(a_matrix: &Matrix, t_var: &Ex) -> bool {
if !a_matrix.is_square() {
return false;
}
let n = a_matrix.nrows();
for i in 0..n {
for j in 0..n {
if a_matrix.get(i, j).contains(t_var) {
return false;
}
}
}
true
}
fn ode_system_is_diagonal(m: &Matrix, n: usize) -> bool {
for i in 0..n {
for j in 0..n {
if i != j && !m.get(i, j).is_zero_structural() {
return false;
}
}
}
true
}
fn solve_ode_system_diagonal(a_matrix: &Matrix, t_var: &Ex, n: usize) -> Vec<Ex> {
let ctx = t_var.context();
(0..n)
.map(|i| {
let ci = ctx.symbol(&format!("C{}", i + 1));
let aii = a_matrix.get(i, i);
if aii.is_zero_structural() {
ci } else {
let exp_term = (aii * t_var).exp();
&ci * &exp_term
}
})
.collect()
}
fn solve_ode_system_expm(a_matrix: &Matrix, t_var: &Ex, n: usize) -> Option<Vec<Ex>> {
let ctx = t_var.context();
let exp_m = a_matrix.matrix_exp_t(t_var).ok()?;
let constants: Vec<Ex> = (1..=n).map(|i| ctx.symbol(&format!("C{i}"))).collect();
let c_vec = Matrix::col_vector(constants).ok()?;
let result = exp_m.matmul(&c_vec).ok()?;
Some((0..n).map(|i| result.get(i, 0).eval()).collect())
}
fn solve_ode_system_eigen(a_matrix: &Matrix, t_var: &Ex, n: usize) -> Option<Vec<Ex>> {
let ctx = t_var.context();
let i_unit = ctx.i_unit();
let eigen = a_matrix.eigenvects().ok()?;
if eigen.iter().map(|(_, _, vs)| vs.len()).sum::<usize>() != n
|| eigen.iter().any(|(_, m, vs)| *m != vs.len())
{
return None;
}
let split = |e: &Ex| -> Option<(Ex, Ex)> {
let (re, im) = e.as_real_imag();
let (re, im) = (re.eval().simplify(), im.eval().simplify());
let clean = |p: &Ex| {
!p.contains(&i_unit) && !p.to_string().contains("re(") && !p.to_string().contains("im(")
};
(clean(&re) && clean(&im)).then_some((re, im))
};
let mut solution: Vec<Ex> = (0..n).map(|_| ctx.int(0)).collect();
let mut const_idx = 1_usize;
for (ev, _, vectors) in &eigen {
let (alpha, beta) = split(ev)?;
if beta.is_zero_structural() {
let exp_ev_t = if alpha.is_zero_structural() {
ctx.int(1)
} else {
(&alpha * t_var).exp()
};
for v in vectors {
let ci = ctx.symbol(&format!("C{const_idx}"));
const_idx += 1;
for (row, sol_row) in solution.iter_mut().enumerate() {
let (vr, vi) = split(v.get(row, 0))?;
if !vi.is_zero_structural() {
return None;
}
if !vr.is_zero_structural() {
*sol_row = &*sol_row + &(&ci * &(&exp_ev_t * &vr));
}
}
}
continue;
}
let b = beta.eval_f64().ok()?;
if b < 0.0 {
continue;
}
let exp_alpha_t = if alpha.is_zero_structural() {
ctx.int(1)
} else {
(&alpha * t_var).exp()
};
let cos_bt = (&beta * t_var).cos();
let sin_bt = (&beta * t_var).sin();
for v in vectors {
let c_a = ctx.symbol(&format!("C{const_idx}"));
let c_b = ctx.symbol(&format!("C{}", const_idx + 1));
const_idx += 2;
for (row, sol_row) in solution.iter_mut().enumerate() {
let (u, w) = split(v.get(row, 0))?;
let m1 = &exp_alpha_t * &(&(&cos_bt * &u) - &(&sin_bt * &w));
let m2 = &exp_alpha_t * &(&(&sin_bt * &u) + &(&cos_bt * &w));
*sol_row = &*sol_row + &(&(&c_a * &m1) + &(&c_b * &m2));
}
}
}
if const_idx != n + 1 {
return None;
}
Some(solution.into_iter().map(|s| s.eval()).collect())
}
#[cfg(test)]
mod tests {
use super::*;
fn sym(a: &mut Arena, name: &str) -> ExprId {
a.symbol(name)
}
fn display(a: &Arena, id: ExprId) -> String {
a.display(id).to_string()
}
#[test]
fn solve_dy_dx_eq_x() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let expr = a.sub(dy, x);
let result = dsolve(&mut a, expr, y, x).expect("should solve");
let s = display(&a, result.solution);
assert!(s.contains("C1"), "should have constant: {s}");
assert!(s.contains("x"), "should contain x: {s}");
}
#[test]
fn solve_dy_dx_eq_0() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let result = dsolve(&mut a, dy, y, x).expect("should solve");
let s = display(&a, result.solution);
assert!(s.contains("C1"), "should be constant: {s}");
}
#[test]
fn solve_second_order_y_plus_y_eq_0() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let d2y = a.intern(ExprNode::Derivative(dy, x));
let expr = a.add(&[d2y, y]);
let r = dsolve(&mut a, expr, y, x).expect("should solve y'' + y = 0");
let s = display(&a, r.solution);
assert!(
s.contains("C1") && s.contains("C2"),
"should have two constants: {s}"
);
assert!(
s.contains("cos") && s.contains("sin"),
"should use trig form (cos and sin): {s}"
);
}
#[test]
fn solve_y_prime_plus_2y_eq_0() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let two = a.int(2);
let two_y = a.mul(&[two, y]);
let expr = a.add(&[dy, two_y]);
let result = dsolve(&mut a, expr, y, x).expect("should solve");
let s = display(&a, result.solution);
assert!(s.contains("C1"), "should have constant: {s}");
assert!(s.contains("exp"), "should contain exp: {s}");
}
#[test]
fn solve_full_separable_xy() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let xy = a.mul(&[x, y]);
let expr = a.sub(dy, xy); let result = dsolve(&mut a, expr, y, x).expect("should solve y' = xy");
let s = display(&a, result.solution);
assert!(s.contains("C1"), "should have constant: {s}");
assert!(s.contains("exp"), "should contain exp: {s}");
}
#[test]
fn solve_variable_coeff_linear_2xy() {
let mut a = Arena::new();
let x = sym(&mut a, "x");
let y = sym(&mut a, "y");
let dy = a.intern(ExprNode::Derivative(y, x));
let two = a.int(2);
let two_x_y = a.mul(&[two, x, y]);
let expr = a.add(&[dy, two_x_y]); let result = dsolve(&mut a, expr, y, x).expect("should solve y' + 2xy = 0");
let s = display(&a, result.solution);
assert!(s.contains("C1"), "should have constant: {s}");
assert!(s.contains("exp"), "should contain exp: {s}");
}
}