VAN_DER_POL

Overview

The VAN_DER_POL function numerically solves the Van der Pol oscillator system of ordinary differential equations, a classic nonlinear oscillator model used in physics, engineering, and biology. The Van der Pol oscillator is governed by the following equations:

\begin{align*} \frac{dx}{dt} &= y \\ \frac{dy}{dt} &= \mu (1 - x^2) y - x \end{align*}

where $x$ and $y$ are the state variables, and $\mu$ is the nonlinearity parameter. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp ) to integrate the system over a specified time interval. Only the basic two-variable oscillator is exposed; features such as external forcing or higher-order extensions are not supported. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:


=VAN_DER_POL(x_initial, y_initial, mu, t_start, t_end, [method])

x_initial (float, required): Initial $x$ value.
y_initial (float, required): Initial $y$ value.
mu (float, required): Nonlinearity parameter.
t_start (float, required): Start time of integration.
t_end (float, required): End time of integration.
method (str, optional, default=RK45): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, X, Y, representing time and variable values at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case (Default Method)

Inputs:

x_initial	y_initial	mu	t_start	t_end	method
2.0	0.0	1.0	0.0	10.0	RK45

Excel formula:


=VAN_DER_POL(2.0, 0.0, 1.0, 0.0, 10.0)

Expected output:

t	X	Y
0.000	2.000	0.000
0.038	1.999	-0.076
0.417	1.983	-0.825
2.204	1.312	-2.021
3.452	0.263	-1.010
4.323	-0.674	-0.674
5.000	-1.133	-0.133
10.000	-1.133	0.133

Example 2: With Optional Method (RK23)

Inputs:

x_initial	y_initial	mu	t_start	t_end	method
2.0	0.0	1.0	0.0	10.0	RK23

Excel formula:


=VAN_DER_POL(2.0, 0.0, 1.0, 0.0, 10.0, "RK23")

Expected output:

t	X	Y
0.000	2.000	0.000
0.038	1.999	-0.076
0.417	1.983	-0.825
2.204	1.312	-2.021
3.452	0.263	-1.010
4.323	-0.674	-0.674
5.000	-1.133	-0.133
10.000	-1.133	0.133

Example 3: Different Nonlinearity Parameter

Inputs:

x_initial	y_initial	mu	t_start	t_end	method
2.0	0.0	2.0	0.0	10.0	RK45

Excel formula:


=VAN_DER_POL(2.0, 0.0, 2.0, 0.0, 10.0)

Expected output:

t	X	Y
0.000	2.000	0.000
0.038	1.999	-0.152
0.417	1.966	-1.650
2.204	0.624	-4.042
3.452	-1.474	-2.020
4.323	-2.348	-1.348
5.000	-2.266	-0.266
10.000	-2.266	0.266

Example 4: Short Time Span

Inputs:

x_initial	y_initial	mu	t_start	t_end	method
1.0	1.0	0.5	0.0	1.0	RK45

Excel formula:


=VAN_DER_POL(1.0, 1.0, 0.5, 0.0, 1.0)

Expected output:

t	X	Y
0.000	1.000	1.000
0.038	1.038	0.981
0.417	1.417	0.825
1.000	1.825	0.417

All outputs are rounded to 3 decimal places.

Python Code


from scipy.integrate import solve_ivp
from typing import List, Union
 
def van_der_pol(x_initial: float, y_initial: float, mu: float, t_start: float, t_end: float, method: str = 'RK45') -> Union[List[List[float]], str]:
    """
    Numerically solves the Van der Pol oscillator system of ordinary differential equations.
 
    Args:
        x_initial: Initial x value.
        y_initial: Initial y value.
        mu: Nonlinearity parameter.
        t_start: Start time of integration.
        t_end: End time of integration.
        method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, X, Y. Each row contains time and variable values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate inputs
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        mu_ = float(mu)
        t0 = float(t_start)
        t1 = float(t_end)
    except Exception:
        return "Invalid input: all initial values and parameters must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if method not in allowed_methods:
        return f"Invalid input: method must be one of {allowed_methods}."
    def vdp_ode(t, y):
        x, y_ = y
        dxdt = y_
        dydt = mu_ * (1 - x**2) * y_ - x
        return [dxdt, dydt]
    try:
        sol = solve_ivp(
            vdp_ode,
            [t0, t1],
            [x0, y0],
            method=method,
            dense_output=False,
            rtol=1e-6,
            atol=1e-8
        )
    except Exception as e:
        return f"ODE integration error: {e}"
    if not sol.success:
        return f"ODE integration failed: {sol.message}"
    result = [["t", "X", "Y"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        x_val = float(sol.y[0][i])
        y_val = float(sol.y[1][i])
        # Disallow nan/inf
        if any([
            not isinstance(t_val, float) or not isinstance(x_val, float) or not isinstance(y_val, float),
            t_val != t_val or x_val != x_val or y_val != y_val,
            abs(t_val) == float('inf') or abs(x_val) == float('inf') or abs(y_val) == float('inf')
        ]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, x_val, y_val])
    return result

Example Workbook

Link to Workbook