Matlab ode45与Scipy odeint求解微分方程组的结果差异及RuntimeWarning问题排查

Hey there, let's break down what's going on with your code and why you're seeing those weird NaN values, plus how to fix it without relying on cranking up the number of sample points.

Your Problem Recap #

You're porting a simple differential equation system from Matlab's ode45 to Scipy's odeint. When using t_list = np.linspace(0, 30, 31) (1-second intervals), you get a RuntimeWarning about invalid values and tons of NaNs in the output. But Matlab's ode45 handles the same model perfectly. When you switch to a denser sample list (linspace(0,30,100)), it works fine. You're wondering two things:

1. Why is this system so sensitive to sample interval?
2. Did you miss any optional parameters in odeint?

First, let's look at your code for reference:

Python Code (Original)

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)

Matlab Code

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

Why Sample Interval Matters (And Why Matlab Works) #

The key difference here is how ode45 and odeint handle step sizes:

• Matlab's ode45 is an adaptive-step solver: It automatically adjusts how big each calculation step is based on how "fast" the system is changing. When your input drops from 1.0 to 0.75 at t=10, or when x/y/z get close to 0 (since you have square roots), ode45 shrinks the step size to avoid overshooting into negative values (which would break the square root and create NaNs).
• Scipy's odeint (which uses the LSODA solver under the hood) is also adaptive-step, but...: When you pass a sparse t_list (like 31 points over 30 seconds), the solver might take larger steps between your specified points. If those steps are too big, it can miss critical changes in the system, leading to x (or y/z) being calculated as negative. Once you have a negative value under a square root, you get NaNs that propagate through all subsequent calculations.

When you use a denser t_list, you're effectively giving the solver more checkpoints, which indirectly limits how large its adaptive steps can be—keeping it from overshooting into negative territory.

Fixing It With odeint Optional Parameters #

You don't need to rely on denser sample points! You can tweak odeint's parameters to match ode45's behavior more closely:

1. Limit the maximum step size with hmax:
   This tells the solver "never take a step larger than X". Setting hmax=0.5 (for example) keeps steps small enough to avoid overshooting, even with a sparse t_list.

2. Tighten error tolerances with rtol and atol:
   The default rtol=1e-3 and atol=1e-6 are pretty loose. Setting stricter values (like rtol=1e-8, atol=1e-10) makes the solver more conservative about error, which helps prevent numerical instability.

Here's your modified Python code with these fixes:

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]
    sys_input = 1.0 if t <= 10 else 0.75
    a = b = c = 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)
yi = [1.0, 1.0, 1.0]
# Add hmax, rtol, and atol to stabilize the solver
ret = odeint(model, y0=yi, t=t_list, hmax=0.5, rtol=1e-8, atol=1e-10)
X = ret[:, 0]
print(X)

# Plot to verify
plt.plot(t_list, X, label='X(t)')
plt.xlabel('Time')
plt.ylabel('X Value')
plt.legend()
plt.show()

This code will give you the same continuous, valid results as Matlab's ode45 without needing a dense sample list.

Quick Summary #

• The NaNs come from negative values being fed into square roots, which happens when odeint takes overly large adaptive steps with a sparse t_list.
• Matlab's ode45 avoids this by default with more aggressive adaptive step sizing.
• You can fix odeint by setting hmax to limit step size, and tightening rtol/atol for better error control.

内容的提问来源于stack exchange，提问作者Ian