Skip to Content

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:

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

where xx and yy 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, [timesteps], [solve_ivp_method])
  • x_initial (float, required): Initial xx value.
  • y_initial (float, required): Initial yy value.
  • mu (float, required): Nonlinearity parameter.
  • t_start (float, required): Start time of integration.
  • t_end (float, required): End time of integration.
  • timesteps (int, optional, default=10): Number of timesteps to solve for.
  • solve_ivp_method (string (enum), optional, default=‘RK45’): Integration method. Valid options: 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_initialy_initialmut_startt_endtimesteps
2.00.01.00.010.010

Excel formula:

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

Output:

tXY
0.0002.0000.000
1.1111.418-0.840
2.222-0.136-2.304
3.333-2.030-0.013
4.444-1.4590.819
5.5560.0312.206
6.6672.0070.065
7.7781.452-0.819
8.889-0.046-2.218
10.000-2.007-0.051

Example 2: With Optional Method (RK23)

Inputs:

x_initialy_initialmut_startt_endtimestepssolve_ivp_method
2.00.01.00.010.010RK23

Excel formula:

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

Output:

tXY
0.0002.0000.000
1.1111.418-0.840
2.222-0.134-2.302
3.333-2.0070.025
4.444-1.4170.843
5.5560.1392.307
6.6672.007-0.029
7.7781.415-0.844
8.889-0.145-2.313
10.000-2.0070.034

Example 3: Different Nonlinearity Parameter

Inputs:

x_initialy_initialmut_startt_endtimesteps
2.00.02.00.010.010

Excel formula:

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

Output:

tXY
0.0002.0000.000
1.1111.652-0.427
2.2220.988-0.908
3.333-1.580-2.796
4.444-1.8390.358
5.556-1.3370.595
6.667-0.0252.569
7.7781.988-0.242
8.8891.581-0.458
10.0000.829-1.106

Example 4: Short Time Span

Inputs:

x_initialy_initialmut_startt_endtimesteps
1.01.00.50.01.010

Excel formula:

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

Output:

tXY
0.0001.0001.000
0.1111.1040.877
0.2221.1940.735
0.3331.2670.580
0.4441.3230.417
0.5561.3600.251
0.6671.3790.088
0.7781.380-0.068
0.8891.364-0.214
1.0001.333-0.352

Python Code

import scipy.integrate as scipy_integrate
import numpy as np
 
def van_der_pol(x_initial, y_initial, mu, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the Van der Pol oscillator system of ordinary differential equations.
 
    This function wraps scipy.integrate.solve_ivp. See:
    https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
 
    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.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_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 input types
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        mu_ = float(mu)
        t0 = float(t_start)
        t1 = float(t_end)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: solve_ivp_method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)
    
    # Van der Pol ODE system
    def vdp_ode(t, y):
        x, y_ = y
        dxdt = y_
        dydt = mu_ * (1 - x**2) * y_ - x
        return [dxdt, dydt]
    # Integrate
    try:
        sol = scipy_integrate.solve_ivp(
            vdp_ode,
            [t0, t1],
            [x0, y0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Integration error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"
    # Format output: header row then data rows
    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,  # NaN check
            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 

Last updated on
SIRSOLVE_BVP