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.

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 (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_evalmethod
-0.50.010.02.00.0RK45
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}, "RK45")

Expected output:

ty
0.02.000
2.00.735
4.00.271
6.00.100
8.00.037
10.00.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.010.00020.000
1.010.17618.665
2.010.16617.589
3.010.03816.689
4.09.83415.916
5.09.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.02.000
2.00.736
4.00.271
6.00.099
8.00.037
10.00.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.010.0005.000
5.05226.2291225.214
10.02927260.501686253.964
15.01638040576.002384014964.530

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 

Last updated on
QUADTRAPEZOID