SOLVE_IVP

Overview

The SOLVE_IVP function numerically integrates a system of ordinary differential equations (ODEs) given an initial value. Specifically, it solves initial value problems (IVPs) of the form:

\frac{dy}{dt} = A \cdot y, \quad y(t_0) = y_0

where A is a coefficient matrix, y is the state vector, and y_0 represents the initial conditions at time t_0. This function wraps SciPy’s solve_ivp, a versatile ODE solver from the scipy.integrate module.

The function supports six integration methods, each suited to different problem characteristics:

RK45 (default): An explicit Runge-Kutta method of order 5(4), based on the Dormand-Prince formula. Recommended for most non-stiff problems.
RK23: A lower-order Runge-Kutta method of order 3(2), using the Bogacki-Shampine formula.
DOP853: A high-precision explicit Runge-Kutta method of order 8, suitable when very accurate solutions are required.
Radau: An implicit Runge-Kutta method of the Radau IIA family (order 5), designed for stiff problems.
BDF: A backward differentiation formula method (variable order 1-5), another implicit method for stiff equations.
LSODA: An Adams/BDF method with automatic stiffness detection and switching, wrapping the classic ODEPACK Fortran solver.

For non-stiff problems, explicit methods (RK45, RK23, DOP853) are preferred. For stiff problems—where the solution contains components varying on vastly different timescales—implicit methods (Radau, BDF, LSODA) provide better stability and efficiency. The LSODA method is particularly useful when stiffness is uncertain, as it automatically detects and adapts to problem characteristics.

All methods use adaptive step-size control, adjusting the integration step to maintain accuracy while minimizing computational cost. The solution is evaluated at user-specified time points within the integration interval.

For more details on Runge-Kutta methods, see Runge-Kutta methods on Wikipedia. The underlying SciPy implementation is available on GitHub.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=SOLVE_IVP(fun_coeffs, t_start, t_end, y_zero, t_evaluate, ivp_method)

fun_coeffs (list[list], required): Coefficient matrix A in the ODE system dy/dt = A·y. Must be square.
t_start (float, required): Initial time for the integration interval.
t_end (float, required): Final time for the integration interval (must be greater than t_start).
y_zero (list[list], required): Initial state vector y(t_start) as a column vector.
t_evaluate (list[list], required): Column vector of times where the solution is computed (must be within [t_start, t_end]).
ivp_method (str, optional, default: “RK45”): Integration method to use.

Returns (list[list]): 2D list [t, y1, y2, …], or error message string.

Examples

Example 1: Simple exponential decay with RK45

Inputs:

fun_coeffs	t_start	t_end	y_zero	t_evaluate
-0.5	0	10	2	0
				2
				4
				6
				8
				10

Excel formula:

=SOLVE_IVP(-0.5, 0, 10, 2, {0;2;4;6;8;10})

Expected output:

Result
0	2
2	0.7353401345874133
4	0.2706911955113017
6	0.09965104748458352
8	0.03669163129975897
10	0.013507816271554396

Example 2: 2D coupled linear system with RK45

Inputs:

fun_coeffs		t_start	t_end	y_zero	t_evaluate	ivp_method
-0.25	0.14	0	5	10	0	RK45
0.25	-0.2			20	1
					2
					3
					4
					5

Excel formula:

=SOLVE_IVP({-0.25,0.14;0.25,-0.2}, 0, 5, {10;20}, {0;1;2;3;4;5}, "RK45")

Expected output:

Result
0	10	20
1	10.1763622202596	18.665189024490083
2	10.165611261041196	17.589409440565607
3	10.038404019834816	16.68858915858497
4	9.83371803690175	15.915760230488328
5	9.582885745113462	15.232945026573965

Example 3: Exponential decay with BDF method

Inputs:

fun_coeffs	t_start	t_end	y_zero	t_evaluate	ivp_method
-0.5	0	10	2	0	BDF
				2
				4
				6
				8
				10

Excel formula:

=SOLVE_IVP(-0.5, 0, 10, 2, {0;2;4;6;8;10}, "BDF")

Expected output:

Result
0	2
2	0.7362429713902789
4	0.27066357188193224
6	0.09947159145323485
8	0.03655917134349607
10	0.013439909120809711

Example 4: 2D system with rapid exponential growth

Inputs:

fun_coeffs		t_start	t_end	y_zero	t_evaluate	ivp_method
1.5	-1	0	15	10	0	RK45
1	-3			5	5
					10
					15

Excel formula:

=SOLVE_IVP({1.5,-1;1,-3}, 0, 15, {10;5}, {0;5;10;15}, "RK45")

Expected output:

Result
0	10	5
5	5226.229049711845	1225.2138594257312
10	2927260.501072503	686253.9635548706
15	1638040576.0015073	384014964.52950907

Python Code

import math
import numpy as np
from scipy.integrate import solve_ivp as scipy_solve_ivp

def solve_ivp(fun_coeffs, t_start, t_end, y_zero, t_evaluate, ivp_method='RK45'):
    """
    Solve an initial value problem for a system of ODEs of the form dy/dt = A @ y.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        fun_coeffs (list[list]): Coefficient matrix A in the ODE system dy/dt = A·y. Must be square.
        t_start (float): Initial time for the integration interval.
        t_end (float): Final time for the integration interval (must be greater than t_start).
        y_zero (list[list]): Initial state vector y(t_start) as a column vector.
        t_evaluate (list[list]): Column vector of times where the solution is computed (must be within [t_start, t_end]).
        ivp_method (str, optional): Integration method to use. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list [t, y1, y2, ...], or error message string.
    """
    def _parse_matrix(data, name):
        """Convert input to matrix format with validation."""
        if isinstance(data, (int, float)):
            if not math.isfinite(data):
                return f"Invalid input: {name} contains non-finite values."
            return [[float(data)]]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        matrix = []
        row_length = None
        for i, row in enumerate(data):
            if not isinstance(row, list) or len(row) == 0:
                return f"Invalid input: row {i+1} in {name} must be a non-empty list."
            try:
                float_row = [float(value) for value in row]
            except (ValueError, TypeError) as e:
                return f"Invalid input: {name} row {i+1} contains non-numeric value."
            if not all(math.isfinite(val) for val in float_row):
                return f"Invalid input: {name} row {i+1} contains non-finite values."
            if row_length is None:
                row_length = len(float_row)
            elif len(float_row) != row_length:
                return f"Invalid input: row {i+1} in {name} has {len(float_row)} columns, expected {row_length}."
            matrix.append(float_row)
        return matrix

    def _parse_column_vector(data, name):
        """Convert input to column vector format with validation."""
        if isinstance(data, (int, float)):
            if not math.isfinite(data):
                return f"Invalid input: {name} contains non-finite value."
            return [float(data)]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        column = []
        for i, row in enumerate(data):
            if not isinstance(row, list) or len(row) != 1:
                return f"Invalid input: row {i+1} in {name} must contain exactly one value."
            try:
                val = float(row[0])
            except (ValueError, TypeError):
                return f"Invalid input: row {i+1} in {name} contains non-numeric value."
            if not math.isfinite(val):
                return f"Invalid input: row {i+1} in {name} contains non-finite value."
            column.append(val)
        return column

    # Validate scalar parameters first (faster checks)
    try:
        start = float(t_start)
        end = float(t_end)
    except (ValueError, TypeError):
        return "Invalid input: t_start and t_end must be numeric values."
    if not math.isfinite(start) or not math.isfinite(end):
        return "Invalid input: t_start and t_end must be finite values."
    if start >= end:
        return "Invalid input: t_start must be strictly less than t_end."

    # Parse and validate coefficient matrix
    matrix_result = _parse_matrix(fun_coeffs, "fun_coeffs")
    if isinstance(matrix_result, str):
        return matrix_result
    coefficient_matrix = np.array(matrix_result, dtype=float)
    if coefficient_matrix.ndim != 2:
        return "Invalid input: fun_coeffs must be a 2D matrix."
    rows, cols = coefficient_matrix.shape
    if rows != cols:
        return f"Invalid input: fun_coeffs must be a square matrix, got {rows}x{cols}."

    # Parse and validate initial state vector
    y0_result = _parse_column_vector(y_zero, "y_zero")
    if isinstance(y0_result, str):
        return y0_result
    initial_state = np.array(y0_result, dtype=float)
    if rows != initial_state.size:
        return f"Invalid input: fun_coeffs dimension ({rows}x{rows}) must match y_zero length ({initial_state.size})."

    # Parse and validate evaluation times
    t_eval_result = _parse_column_vector(t_evaluate, "t_evaluate")
    if isinstance(t_eval_result, str):
        return t_eval_result
    if len(t_eval_result) == 0:
        return "Invalid input: t_evaluate must contain at least one time point."
    evaluation_times = np.array(t_eval_result, dtype=float)

    # Validate evaluation times ordering and bounds
    if len(evaluation_times) > 1 and not np.all(np.diff(evaluation_times) >= 0):
        return "Invalid input: t_evaluate must be sorted in non-decreasing order."
    if evaluation_times[0] < start:
        return f"Invalid input: first t_evaluate point ({evaluation_times[0]}) is before t_start ({start})."
    if evaluation_times[-1] > end:
        return f"Invalid input: last t_evaluate point ({evaluation_times[-1]}) is after t_end ({end})."

    # Validate and normalize method parameter
    try:
        method_key = str(ivp_method).strip().upper()
    except (ValueError, TypeError, AttributeError):
        return "Invalid input: ivp_method must be a string."

    allowed_methods = {
        "RK45": "RK45",
        "RK23": "RK23",
        "DOP853": "DOP853",
        "RADAU": "Radau",
        "BDF": "BDF",
        "LSODA": "LSODA",
    }
    if method_key not in allowed_methods:
        valid_methods = ", ".join(["RK45", "RK23", "DOP853", "Radau", "BDF", "LSODA"])
        return f"Invalid input: ivp_method must be one of: {valid_methods} (got '{ivp_method}')."
    method_value = allowed_methods[method_key]

    # Define the ODE system
    def ode_system(_t: float, state: np.ndarray) -> np.ndarray:
        """Linear system dy/dt = A·y."""
        return coefficient_matrix @ state

    # Solve the initial value problem
    try:
        solution = scipy_solve_ivp(
            ode_system,
            (start, end),
            initial_state,
            method=method_value,
            t_eval=evaluation_times,
        )
    except Exception as exc:  # noqa: BLE001
        return f"Integration failed: {type(exc).__name__}: {exc}"

    if not solution.success:
        return f"Integration failed: {solution.message}"

    # Combine time and solution into output format
    combined = np.hstack((solution.t[:, np.newaxis], solution.y.T))
    if not np.all(np.isfinite(combined)):
        return "Integration failed: result contains non-finite values (NaN or Inf)."

    return combined.tolist()

Overview

Excel Usage

Examples

Python Code

Online Calculator