Differences between Matlab ode45 and Scipy odeint: same model different results

Question

EDIT

So I was able to produce the correct plot, but it was only after adjusting the sampling interval to the following:

t_list = np.linspace(0, 30, 100)

which prints X as:

[1.         1.         1.         1.         1.         1.
 1.         1.         1.         1.         1.         1.
 1.         1.         1.         1.         1.         1.
 0.91265299 0.8107059  0.7366542  0.68370578 0.64633005 0.62021062
 0.6020953  0.58960093 0.58101775 ...]

But that begs the question: why is this system so dependent on the sampling intervals?

END EDIT

I am trying to recreate a simple matlab system of differential equation into python using scipy. I do not know why I am receiving a RuntimeWarning: invalid value encountered in double_scalars in the scipy execution. Am I missing an optional parameter in the odeint call?

python:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def model(y0, t):
    x = y0[0]
    y = y0[1]
    z = y0[2]
    if t <= 10:
        sys_input = 1.0
    else:
        sys_input = 0.75
    a = 1.0
    b = 1.0
    c = 1.0
    E = 1.0
    dxdt = sys_input - a * E * (x ** 0.5)
    dydt = a * E * (x ** 0.5) - b * (y ** 0.5)
    dzdt = b * (y ** 0.5) - c * (z ** 0.5)
    return [dxdt, dydt, dzdt]


t_list = np.linspace(0, 30, 31)

# Initial conditions vector
yi = [1.0, 1.0, 1.0]
ret = odeint(model, y0=yi, t=t_list)
X = ret[:, 0]
print(X)

and it prints the following:

<input>:18: RuntimeWarning: invalid value encountered in double_scalars
<input>:19: RuntimeWarning: invalid value encountered in double_scalars
[ 1.  1.  1.  1. nan nan nan nan nan nan nan nan nan nan nan nan nan nan
 nan nan nan nan nan nan nan nan nan nan nan nan nan]

where as the following matlab code produces a continuous result:

tspan = [0,30];
x0 = 1.0; % Initial value of x
y0 = 1.0; % Initial value of y
z0 = 1.0; % Initial value of z
initial_values = [x0; y0; z0]; % Initial value of the vector w
[T,R] = ode45(@(t,w) diff_eq(t,w),tspan,initial_values);
X = R(:,1);
Y = R(:,2);
Z = R(:,3);

for i = 1: length(X)
if(mod(i, 10)==0 &&  i > 1)
    disp(' ');
end
fprintf('X[%i] = %.2f, ', i, X(i));
end
disp(' ');


function dw_vectordt = diff_eq(t,w_vector) 
 x = w_vector(1);
 y = w_vector(2);
 z = w_vector(3);
 if (t<=10)
    sys_input= 1.0;    
else
    sys_input=0.75;
end
 a = 1.0;
 b = 1.0;
 c = 1.0;
 E = 1.0;
 dxdt = sys_input-a*E*x^(0.5);
 dydt = a*E*x^(0.5)-b*y^(0.5);
 dzdt = b*y^(0.5)-c*z^(0.5);
 dw_vectordt = [dxdt; dydt; dzdt];
end

the print statement returns:

X[1] = 1.00, X[2] = 1.00, X[3] = 1.00, X[4] = 1.00, X[5] = 1.00, X[6] = 1.00, X[7] = 1.00, X[8] = 1.00, X[9] = 1.00,  
X[10] = 1.00, X[11] = 1.00, X[12] = 1.00, X[13] = 1.00, X[14] = 1.01, X[15] = 0.99, X[16] = 0.92, X[17] = 0.82, X[18] = 0.78, X[19] = 0.74,  
X[20] = 0.71, X[21] = 0.68, X[22] = 0.66, X[23] = 0.64, X[24] = 0.63, X[25] = 0.61, X[26] = 0.60, X[27] = 0.59, X[28] = 0.59, X[29] = 0.58,  
X[30] = 0.58, X[31] = 0.57, X[32] = 0.57, X[33] = 0.57, X[34] = 0.57, X[35] = 0.57, X[36] = 0.56, X[37] = 0.56, X[38] = 0.56, X[39] = 0.56,  
X[40] = 0.56, X[41] = 0.56, X[42] = 0.56, X[43] = 0.56, X[44] = 0.56, X[45] = 0.56, X[46] = 0.56, X[47] = 0.56, X[48] = 0.56, X[49] = 0.56,  
X[50] = 0.56, X[51] = 0.56, X[52] = 0.56, X[53] = 0.56,

Answer 1

The ODE system you are dealing with is likely stiff. The RuntimeWarning you are encountering is raised by the square root operations as elements of y0 become negative during integration. This is occurring because the timestep of the integrator is too large, and the solver you are using is not well-suited for potentially stiff systems. Increasing the number of elements provided to t through t_list likely decreases the timestep and, therefore, allows for a workaround. To better understand what is happening, I encourage you to play with the snippet below, which uses the newer solve_ivp API, recommended by SciPy. Of particular interest are the keyword arguments method and max_step. Using either of RK23, RK45, DOP853 and LSODA as a method will result in poor solution estimates out of the box as all of these solvers start with a non-stiff method. LSODA should detect the stiffness and switch to a stiff integrator, but it fails to do so as the timestep quickly increases. For all of these methods, setting max_step=0.5 allows dealing with the ODE system at the potential cost of computation time. Alternatively, using either Radau or BDF will work out of the box as these solvers can deal with stiff ODE systems. However, it is recommended to manually provide the Jacobian of the system, otherwise it is approximated by finite differences.

import numpy as np
from scipy.integrate import solve_ivp


def model(t, y0):
    x, y, z = y0
    sys_input = 1 if t <= 10 else 0.75
    a, b, c, E = 1, 1, 1, 1
    return (sys_input - a * E * np.sqrt(x),
            a * E * np.sqrt(x) - b * np.sqrt(y),
            b * np.sqrt(y) - c * np.sqrt(z))


sol = solve_ivp(model, t_span=(0, 30), y0=(1, 1, 1))
print(sol)