mathcore::differential — ODE and PDE Solvers

The mathcore::differential module provides numerical solvers for ordinary and partial differential equations. DifferentialEquations handles scalar first-order ODEs (RK4), systems of ODEs, second-order ODEs, Euler’s method, stiff ODEs (implicit Euler with Newton iteration), and linear constant-coefficient ODEs. The companion PDESolver struct implements finite-difference schemes for the heat equation, wave equation, and Laplace equation. All time-evolution solvers return an ODESolution containing sampled time and state values.

`pub struct ODESolution`

The result type returned by all ODE solver methods.

#[derive(Debug, Clone)]
pub struct ODESolution {
    pub t: Vec<f64>,
    pub y: Vec<Vec<f64>>,
}

Field	Type	Description
`t`	`Vec<f64>`	Monotonically increasing time (or spatial) values, starting at `t0` and ending at `t_final`. Length is `steps + 1`.
`y`	`Vec<Vec<f64>>`	State values at each time step. Each inner `Vec<f64>` holds the values of all dependent variables. For scalar ODEs the inner vector has length 1; for systems it has length equal to the number of equations.

`pub enum BoundaryCondition`

Encodes the two standard forms of boundary data for ODE problems.

#[derive(Debug, Clone)]
pub enum BoundaryCondition {
    /// Initial value problem: state y0 known at time t0
    InitialValue {
        t0: f64,
        y0: Vec<f64>,
    },
    /// Boundary value problem: partial state known at both endpoints
    BoundaryValue {
        ta: f64,
        tb: f64,
        ya: Vec<f64>,
        yb: Vec<f64>,
    },
}

Variant	Fields	Use case
`InitialValue`	`t0`, `y0`	Standard Cauchy / IVP problems where all state is known at the start.
`BoundaryValue`	`ta`, `tb`, `ya`, `yb`	BVPs where partial state is specified at both ends of the domain.

`pub struct DifferentialEquations`

A zero-sized unit struct grouping all ODE solver methods.

pub struct DifferentialEquations;

`solve_ode_first_order`

Solves a scalar first-order ODE dy/dt = f(t, y) using the fourth-order Runge-Kutta (RK4) method. The expression f is evaluated symbolically at each stage of every step by substituting the current t and y values via the engine. RK4 update rule (per step of size h):

k1 = f(t,       y)
k2 = f(t + h/2, y + h*k1/2)
k3 = f(t + h/2, y + h*k2/2)
k4 = f(t + h,   y + h*k3)
y  ← y + h*(k1 + 2*k2 + 2*k3 + k4)/6

pub fn solve_ode_first_order(
    f: &Expr,
    t_var: &str,
    y_var: &str,
    initial_condition: (f64, f64),
    t_final: f64,
    steps: usize,
) -> Result<ODESolution, MathError>

&Expr

required

The right-hand side expression f(t, y). Must evaluate to Expr::Number when t_var and y_var are substituted with numeric values.

t_var

&str

required

The name of the independent variable (time). Bound at each evaluation.

y_var

&str

required

The name of the dependent variable (state). Bound at each stage of RK4.

initial_condition

(f64, f64)

required

Tuple (t0, y0) — the initial time and initial state value.

t_final

f64

required

The end time of integration. The step size is computed as h = (t_final - t0) / steps.

steps

usize

required

Number of equal sub-intervals. The solution grid has steps + 1 points.

Returns Ok(ODESolution) with t and y vectors, or Err(MathError::InvalidOperation) if f evaluates to a non-numeric result at any stage.

use mathcore::differential::DifferentialEquations;
use mathcore::parser::Parser;

// dy/dt = -y,  y(0) = 1  →  y(t) = e^{-t}
let f = Parser::parse("-y").unwrap();
let sol = DifferentialEquations::solve_ode_first_order(
    &f, "t", "y",
    (0.0, 1.0), // y(0) = 1
    1.0,        // integrate to t = 1
    100,        // 100 steps
).unwrap();

let y_at_1 = sol.y.last().unwrap()[0];
println!("y(1) ≈ {:.6}", y_at_1); // ≈ 0.367879  (e^{-1})
assert!((y_at_1 - 0.3679).abs() < 0.01);

`euler_method`

Solves dy/dt = f(t, y) using the simple first-order forward Euler method. Less accurate than RK4 but computationally cheaper and useful for pedagogical comparison.

pub fn euler_method(
    f: &Expr,
    t_var: &str,
    y_var: &str,
    initial_condition: (f64, f64),
    t_final: f64,
    steps: usize,
) -> Result<ODESolution, MathError>

Parameters and return value are identical to solve_ode_first_order. The update rule is y ← y + h * f(t, y).

`solve_stiff_ode`

Solves a stiff first-order scalar ODE using the implicit (backward) Euler method with Newton iteration to resolve the implicit equation at each step. The jacobian parameter optionally supplies a symbolic expression for ∂f/∂y; when None, a numerical finite-difference approximation is used.

pub fn solve_stiff_ode(
    f: &Expr,
    jacobian: Option<&Expr>,
    t_var: &str,
    y_var: &str,
    initial_condition: (f64, f64),
    t_final: f64,
    steps: usize,
    tolerance: f64,
) -> Result<ODESolution, MathError>

&Expr

required

The right-hand side expression f(t, y).

jacobian

Option<&Expr>

required

Optional expression for ∂f/∂y. When None, the Jacobian is approximated numerically with a finite-difference step of 1e-8.

t_var

&str

required

Independent variable name.

y_var

&str

required

Dependent variable name.

initial_condition

(f64, f64)

required

Tuple (t0, y0).

t_final

f64

required

End time.

steps

usize

required

Number of time steps.

tolerance

f64

required

Newton-iteration convergence tolerance. The inner Newton loop runs up to 10 iterations; it terminates early when |residual| < tolerance.

`solve_ode_second_order`

Reduces a second-order linear ODE y'' + p(x)*y' + q(x)*y = r(x) to a system of two first-order equations and solves it with RK4.

pub fn solve_ode_second_order(
    p: &Expr,
    q: &Expr,
    r: &Expr,
    x_var: &str,
    initial_conditions: (f64, f64, f64),
    x_final: f64,
    steps: usize,
) -> Result<ODESolution, MathError>

&Expr

required

Coefficient of y' (may depend on x_var).

&Expr

required

Coefficient of y (may depend on x_var).

&Expr

required

The forcing (right-hand side) function (may depend on x_var).

x_var

&str

required

The name of the independent variable.

initial_conditions

(f64, f64, f64)

required

Tuple (x0, y(x0), y'(x0)) — starting position, value, and first derivative.

x_final

f64

required

End of the integration domain.

steps

usize

required

Number of steps. Each y vector in the solution contains [y, y'] at the corresponding x.

`solve_ode_system`

Solves a system of n coupled first-order ODEs dy_i/dt = f_i(t, y₁, ..., y_n) simultaneously using RK4.

pub fn solve_ode_system(
    functions: &[Expr],
    t_var: &str,
    y_vars: &[String],
    initial_conditions: (f64, Vec<f64>),
    t_final: f64,
    steps: usize,
) -> Result<ODESolution, MathError>

functions

&[Expr]

required

Slice of expressions f₁, ..., fₙ. Must have the same length as y_vars and initial_conditions.1.

t_var

&str

required

The independent variable name.

y_vars

&[String]

required

Names of the dependent variables in the same order as functions.

initial_conditions

(f64, Vec<f64>)

required

Tuple (t0, [y₁(t0), y₂(t0), ..., yₙ(t0)]).

t_final

f64

required

Integration end time.

steps

usize

required

Number of time steps. Each y inner vector has length n.

use mathcore::differential::DifferentialEquations;
use mathcore::parser::Parser;

// Circular motion:  dx/dt = -y,  dy/dt = x
// x(0) = 1, y(0) = 0  →  x ≈ cos(t), y ≈ sin(t)
let f1 = Parser::parse("-y").unwrap();
let f2 = Parser::parse("x").unwrap();

let sol = DifferentialEquations::solve_ode_system(
    &[f1, f2],
    "t",
    &["x".to_string(), "y".to_string()],
    (0.0, vec![1.0, 0.0]),
    6.283185, // 2π
    1000,
).unwrap();

let final_x = sol.y.last().unwrap()[0];
let final_y = sol.y.last().unwrap()[1];
println!("x(2π) ≈ {:.4}", final_x); // ≈ 1.0
println!("y(2π) ≈ {:.4}", final_y); // ≈ 0.0

`solve_linear_constant_coeff`

Finds the general analytical solution of a linear ODE with constant coefficients by factoring the characteristic polynomial and constructing the basis of exponential (and oscillatory) solutions. The equation is represented by its coefficient vector [aₙ, aₙ₋₁, ..., a₁, a₀] for the operator aₙ y^(n) + ... + a₁ y' + a₀ y = 0. The method builds the characteristic polynomial aₙ rⁿ + ... + a₁ r + a₀, finds its roots via Solver::solve, and assembles the general solution using terms of the form Cᵢ · e^(rᵢ t) for real roots and e^(αt) · (Cᵢ cos(βt) + Cⱼ sin(βt)) for complex-conjugate pairs.

pub fn solve_linear_constant_coeff(
    coefficients: &[f64],
    initial_conditions: &[f64],
) -> Result<Expr, MathError>

coefficients

&[f64]

required

Coefficients of the ODE in descending order of derivative: [aₙ, aₙ₋₁, ..., a₁, a₀]. The leading coefficient a₀ corresponds to the highest-order derivative y^(n). Must be non-empty.

initial_conditions

&[f64]

required

Initial values [y(0), y'(0), ..., y^(n-1)(0)]. The length must equal coefficients.len() - 1 (i.e. the order of the ODE).

Returns Ok(Expr) containing the symbolic general solution (with undetermined constants C1, C2, etc.), or Err(MathError) if the coefficient or initial condition arrays are empty or have mismatched sizes.

use mathcore::differential::DifferentialEquations;

// y'' + y = 0  (coefficients [1, 0, 1], order 2)
// General solution: C1*cos(t) + C2*sin(t)
let sol = DifferentialEquations::solve_linear_constant_coeff(
    &[1.0, 0.0, 1.0], // a2=1, a1=0, a0=1
    &[1.0, 0.0],      // y(0) = 1, y'(0) = 0
).unwrap();
println!("{}", sol);

`pub struct PDESolver`

Implements finite-difference PDE solvers for the three classical second-order PDEs.

pub struct PDESolver;

`solve_heat_equation`

Solves the 1-D heat equation ∂u/∂t = α ∂²u/∂x² using the explicit (forward-in-time, centred-in-space) finite-difference scheme. Enforces the CFL-like stability criterion α * dt/dx² ≤ 0.5; violating it returns an error.

pub fn solve_heat_equation(
    alpha: f64,
    initial_condition: &dyn Fn(f64) -> f64,
    boundary_left: f64,
    boundary_right: f64,
    x_range: (f64, f64),
    t_final: f64,
    nx: usize,
    nt: usize,
) -> Result<Vec<Vec<f64>>, MathError>

alpha

f64

required

Thermal diffusivity coefficient.

initial_condition

&dyn Fn(f64) -> f64

required

Closure evaluated at each spatial grid point to set u(x, 0).

boundary_left

f64

required

Dirichlet boundary value at x = x_min for all t > 0.

boundary_right

f64

required

Dirichlet boundary value at x = x_max for all t > 0.

x_range

(f64, f64)

required

Spatial domain (x_min, x_max).

t_final

f64

required

Total simulation time.

usize

required

Number of spatial grid points.

usize

required

Number of time steps.

Returns Ok(Vec<Vec<f64>>) of shape nt × nx containing the field u(t, x), or Err if the stability criterion is violated.

`solve_wave_equation`

Solves the 1-D wave equation ∂²u/∂t² = c² ∂²u/∂x² using an explicit second-order finite-difference scheme. Enforces the CFL condition c * dt/dx ≤ 1.

pub fn solve_wave_equation(
    c: f64,
    initial_position: &dyn Fn(f64) -> f64,
    initial_velocity: &dyn Fn(f64) -> f64,
    x_range: (f64, f64),
    t_final: f64,
    nx: usize,
    nt: usize,
) -> Result<Vec<Vec<f64>>, MathError>

f64

required

Wave propagation speed.

initial_position

&dyn Fn(f64) -> f64

required

Sets u(x, 0) — the initial displacement profile.

initial_velocity

&dyn Fn(f64) -> f64

required

Sets ∂u/∂t(x, 0) — the initial velocity profile.

`solve_laplace_equation`

Solves the 2-D Laplace equation ∇²u = 0 on a rectangular domain using Gauss-Seidel iteration. Boundary values are supplied via a closure that returns Some(value) at boundary points and None at interior points.

pub fn solve_laplace_equation(
    boundary_conditions: &dyn Fn(f64, f64) -> Option<f64>,
    x_range: (f64, f64),
    y_range: (f64, f64),
    nx: usize,
    ny: usize,
    tolerance: f64,
    max_iterations: usize,
) -> Result<Vec<Vec<f64>>, MathError>

boundary_conditions

&dyn Fn(f64, f64) -> Option<f64>

required

Given coordinates (x, y), returns Some(u) to fix the value as a boundary condition, or None to mark the point as an interior unknown.

tolerance

f64

required

Gauss-Seidel convergence threshold. Iteration stops early when the maximum change in any cell drops below this value.

max_iterations

usize

required

Hard upper bound on iteration count regardless of convergence.

Core API

Modules

mathcore::differential — ODE and PDE Solvers

`pub struct ODESolution`

`pub enum BoundaryCondition`

`pub struct DifferentialEquations`

`solve_ode_first_order`

`euler_method`

`solve_stiff_ode`

`solve_ode_second_order`

`solve_ode_system`

`solve_linear_constant_coeff`

`pub struct PDESolver`

`solve_heat_equation`

`solve_wave_equation`

`solve_laplace_equation`

Build docs developers (and LLMs) love

Core API

Modules

Documentation Index

​pub struct ODESolution

​pub enum BoundaryCondition

​pub struct DifferentialEquations

​solve_ode_first_order

​euler_method

​solve_stiff_ode

​solve_ode_second_order

​solve_ode_system

​solve_linear_constant_coeff

​pub struct PDESolver

​solve_heat_equation

​solve_wave_equation

​solve_laplace_equation

Build docs developers (and LLMs) love

`pub struct ODESolution`

`pub enum BoundaryCondition`

`pub struct DifferentialEquations`

`solve_ode_first_order`

`euler_method`

`solve_stiff_ode`

`solve_ode_second_order`

`solve_ode_system`

`solve_linear_constant_coeff`

`pub struct PDESolver`

`solve_heat_equation`

`solve_wave_equation`

`solve_laplace_equation`