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:

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

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

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

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): A 2D list representing the matrix A in the ODE system dy/dt = A @ y.
t_start (float, required): The initial time t0 for the integration interval.
t_end (float, required): The final time tf for the integration interval.
y_zero (2D list, required): A 2D list (column vector) representing the initial state y(t0).
t_eval (2D list, required): A 2D list (column vector) of times at which to store the computed solution. Must be sorted and lie within t_start and t_end.
method (string, optional, default=‘RK45’): The integration method to use. Supported methods are ‘RK45’, ‘RK23’, ‘DOP853’, ‘Radau’, ‘BDF’, ‘LSODA’.

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

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

Inputs:

Python Code


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.
 
    Args:
        fun_coeffs: 2D list, representing the matrix A in dy/dt = A @ y.
        t_start: float, the initial time t0 for the integration interval.
        t_end: float, the final time tf for the integration interval.
        y_zero: 2D list (column vector), representing the initial state y(t0).
        t_evaluate: 2D list (column vector), of times at which to store the computed solution.
        method: string, the integration method to use (default: 'RK45').
 
    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 (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Handle Excel scalar case for fun_coeffs
    if isinstance(fun_coeffs, (int, float)):
        fun_coeffs = [[float(fun_coeffs)]]
    # Validate fun_coeffs
    if not isinstance(fun_coeffs, list) or not all(isinstance(row, list) for row in fun_coeffs):
        return "Invalid input: fun_coeffs must be a 2D list."
    try:
        A = np.array(fun_coeffs, dtype=float)
    except ValueError:
        return "Invalid input: fun_coeffs must contain numeric values."
    if A.ndim != 2:
        return "Invalid input: fun_coeffs must be a 2D list (matrix)."
 
    # Validate t_start and t_end
    try:
        t_span = [float(t_start), float(t_end)]
    except ValueError:
        return "Invalid input: t_start and t_end must be numbers."
    if t_span[0] >= t_span[1]:
        return "Invalid input: t_start must be less than t_end."
 
    # Handle Excel scalar case for y_zero
    if isinstance(y_zero, (int, float)):
        y_zero = [[float(y_zero)]]
    # Validate y_zero
    if not isinstance(y_zero, list) or not all(isinstance(row, list) for row in y_zero):
        return "Invalid input: y_zero must be a 2D list (column vector)."
    try:
        y0_np = np.array(y_zero, dtype=float).flatten()
    except ValueError:
        return "Invalid input: y_zero must contain numeric values."
    if y0_np.ndim != 1:
        return "Invalid input: y_zero must be a 2D list representing a column vector."
    if A.shape[0] != A.shape[1] or A.shape[0] != len(y0_np):
        return "Invalid input: fun_coeffs matrix must be square and its dimension must match the length of y0."
 
    # Validate t_evaluate
    if not isinstance(t_evaluate, list) or not all(isinstance(row, list) for row in t_evaluate):
        return "Invalid input: t_evaluate must be a 2D list (column vector)."
    try:
        t_eval_np = np.array(t_evaluate, dtype=float).flatten()
    except ValueError:
        return "Invalid input: t_evaluate must contain numeric values."
    if t_eval_np.ndim != 1:
        return "Invalid input: t_evaluate must be a 2D list representing a column vector."
    if not np.all(np.diff(t_eval_np) >= 0):
        return "Invalid input: t_evaluate must be sorted in non-decreasing order."
    if not (t_span[0] <= t_eval_np[0] and t_eval_np[-1] <= t_span[1]):
        return "Invalid input: t_eval points must lie within t_start and t_end."
 
    # Validate method
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if method not in allowed_methods:
        return f"Invalid input: method must be one of {', '.join(allowed_methods)}."
 
    # Define the ODE function
    def fun(t, y):
        return A @ y
 
    try:
        sol = scipy_solve_ivp(fun, t_span, y0_np, method=method, t_eval=t_eval_np)
        if not sol.success:
            return f"Solver failed: {sol.message}"
 
        # Combine t_eval and y values into a single 2D list
        output = np.hstack((sol.t[:, np.newaxis], sol.y.T))
        return output.tolist()
    except Exception as e:
        return f"An unexpected error occurred during integration: {e}"

Live Demo

Example Workbook

Link to Workbook