Skip to Content

SOLVE_IVP

Overview

The SOLVE_IVP function numerically integrates a system of ordinary differential equations (ODEs) given an initial value. It is a wrapper around the scipy.integrate.solve_ivp method, simplified for use in Excel.

The general form of an initial value problem for a system of ODEs is:

dydt=f(t,y)y(t0)=y0\frac{dy}{dt} = f(t, y) y(t_0) = y_0

Here, tt is a 1-D independent variable (time), y(t)y(t) is an N-D vector-valued function (state), and f(t,y)f(t, y) is an N-D vector-valued function that determines the differential equations. The goal is to find y(t)y(t) approximately satisfying the differential equations, given an initial value y(t0)=y0y(t_0)=y_0.

This function simplifies the f(t, y) callable by assuming a linear system of the form dy/dt=A⋅ydy/dt = A \cdot y, where AA is a matrix provided as fun_coeffs. This means the right-hand side function f does not depend on tt explicitly and is linear in y. Only the most commonly used parameters are exposed, mirroring the defaults from scipy.integrate.solve_ivp.

For more details on the underlying Scipy function, refer to the official SciPy documentation .

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

Usage

To use the function in Excel:

=SOLVE_IVP(fun_coeffs, t_start, t_end, y_zero, t_eval, [method])
  • fun_coeffs (2D list, required): Square matrix A in dy/dt = A @ y with numeric entries.
  • t_start (float, required): Initial time t0 for the integration interval.
  • t_end (float, required): Final time tf for the integration interval; must satisfy t_end > t_start.
  • y_zero (2D list, required): Column vector representing y(t0) with exactly one value per row.
  • t_eval (2D list, required): Column vector of evaluation times containing at least one value, sorted in non-decreasing order, and within [t_start, t_end].
  • method (string, optional, default=‘RK45’): Solver to use; supported values are RK45, RK23, DOP853, Radau, BDF, and LSODA (case-insensitive).

The function returns a 2D list (matrix) where the first column is the time points from t_eval and subsequent columns are the corresponding values of the solution y at those time points. Returns an error message (string) if inputs are invalid.

Examples

Example 1: Simple Exponential Decay (RK45)

This example solves the simple exponential decay ODE dy/dt=−0.5ydy/dt = -0.5y with an initial condition y(0)=2y(0)=2 using the RK45 method. The solution is evaluated at specific time points.

Inputs:

fun_coeffst_startt_endy_zerot_eval
-0.50.010.02.00.0
2.0
4.0
6.0
8.0
10.0

Excel formula:

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

Expected output:

ty
0.0002.000
2.0000.735
4.0000.271
6.0000.100
8.0000.037
10.0000.014

Example 2: 2D Linear System (RK45)

This example solves a 2D system dy/dt=Aydy/dt = A y with A=[−0.250.140.25−0.2]A = \begin{bmatrix}-0.25 & 0.14 \\ 0.25 & -0.2\end{bmatrix} and initial state y(0)=[10,20]Ty(0) = [10, 20]^T using the RK45 method.

Inputs:

fun_coeffst_startt_endy_zerot_evalmethod
-0.250.140.05.010.020.00.0RK45
0.25-0.21.0
2.0
3.0
4.0
5.0

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:

ty1y2
0.00010.00020.000
1.00010.17618.665
2.00010.16617.589
3.00010.03816.689
4.0009.83415.916
5.0009.58315.233

Example 3: Exponential Decay (BDF method)

This example solves the same scalar decay ODE as Example 1, but using the BDF method.

Inputs:

fun_coeffst_startt_endy_zerot_evalmethod
-0.50.010.02.00.0BDF
2.0
4.0
6.0
8.0
10.0

Excel formula:

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

Expected output:

ty
0.0002.000
2.0000.736
4.0000.271
6.0000.099
8.0000.037
10.0000.013

Example 4: 2D System with Large Growth (RK45)

This example solves a 2D system dy/dt=Aydy/dt = A y with A=[1.5−1.01.0−3.0]A = \begin{bmatrix}1.5 & -1.0 \\ 1.0 & -3.0\end{bmatrix} and initial state y(0)=[10,5]Ty(0) = [10, 5]^T over a longer interval, showing rapid growth.

Inputs:

fun_coeffst_startt_endy_zerot_evalmethod
1.5-1.00.015.010.05.00.0RK45
1.0-3.05.0
10.0
15.0

Excel formula:

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

Expected output:

ty1y2
0.00010.0005.000
5.0005226.2291225.214
10.0002927260.501686253.964
15.0001638040576.002384014964.530

Python Code

from typing import List, Sequence
 
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, method="RK45"):
    """
    Solve an initial value problem for a system of ODEs of the form dy/dt = A @ y.
 
    This function is a lightweight wrapper around `scipy.integrate.solve_ivp`.
    For details see: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
 
    Args:
        fun_coeffs (2D list): Matrix A in dy/dt = A @ y. Must be square with numeric entries.
        t_start (float): Initial time of integration.
        t_end (float): Final time of integration (must be greater than t_start).
        y_zero (2D list): Column vector for the initial state y(t_start); must have at least one row.
        t_evaluate (2D list): Column vector of times (within [t_start, t_end]) where the solution is stored.
        method (str, optional): Integration method name ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA').
 
    Returns:
        2D list where the first column contains evaluation times and remaining columns contain solution values,
        or an error message (str) when inputs are invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
 
    def _ensure_matrix(data: Sequence[Sequence[float | int]] | float | int, name: str) -> List[List[float]] | str:
        if isinstance(data, (int, float)):
            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: List[List[float]] = []
        row_length: int | None = None
        for row in data:
            if not isinstance(row, list) or len(row) == 0:
                return f"Invalid input: each row in {name} must be a non-empty list."
            try:
                float_row = [float(value) for value in row]
            except Exception:
                return f"Invalid input: {name} must contain numeric values."
            if row_length is None:
                row_length = len(float_row)
            elif len(float_row) != row_length:
                return f"Invalid input: all rows in {name} must have the same length."
            matrix.append(float_row)
        return matrix
 
    def _ensure_column(data: Sequence[Sequence[float | int]] | float | int, name: str) -> List[float] | str:
        if isinstance(data, (int, float)):
            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: List[float] = []
        for row in data:
            if not isinstance(row, list) or len(row) != 1:
                return f"Invalid input: each row in {name} must contain exactly one value."
            try:
                column.append(float(row[0]))
            except Exception:
                return f"Invalid input: {name} must contain numeric values."
        return column
 
    matrix_result = _ensure_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 list (matrix)."
    if coefficient_matrix.shape[0] != coefficient_matrix.shape[1]:
        return "Invalid input: fun_coeffs must form a square matrix."
 
    y0_result = _ensure_column(y_zero, "y_zero")
    if isinstance(y0_result, str):
        return y0_result
    initial_state = np.array(y0_result, dtype=float)
    if coefficient_matrix.shape[0] != initial_state.size:
        return "Invalid input: fun_coeffs dimension must match the length of y_zero."
 
    t_eval_result = _ensure_column(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)
 
    try:
        start = float(t_start)
        end = float(t_end)
    except Exception:
        return "Invalid input: t_start and t_end must be numbers."
    if not math.isfinite(start) or not math.isfinite(end):
        return "Invalid input: t_start and t_end must be finite."
    if start >= end:
        return "Invalid input: t_start must be less than t_end."
 
    if 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 or evaluation_times[-1] > end:
        return "Invalid input: t_evaluate points must lie within t_start and t_end."
 
    try:
        method_key = str(method).strip().upper()
    except Exception:
        return "Invalid input: 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:
        ordered_methods = ", ".join(["RK45", "RK23", "DOP853", "Radau", "BDF", "LSODA"])
        return f"Invalid input: method must be one of {ordered_methods}."
    method_value = allowed_methods[method_key]
 
    def system(_t: float, state: np.ndarray) -> np.ndarray:
        return coefficient_matrix @ state
 
    try:
        solution = scipy_solve_ivp(
            system,
            (start, end),
            initial_state,
            method=method_value,
            t_eval=evaluation_times,
        )
    except Exception as exc:  # noqa: BLE001
        return f"scipy.integrate.solve_ivp error: {exc}"
 
    if not solution.success:
        return f"scipy.integrate.solve_ivp error: {solution.message}"
 
    combined = np.hstack((solution.t[:, np.newaxis], solution.y.T))
    if not np.all(np.isfinite(combined)):
        return "scipy.integrate.solve_ivp error: result is not finite."
 
    return combined.tolist()

Example Workbook

Link to Workbook 

Last updated on
SOLVE_BVPHESSIAN