Matlab和ode45的差异:相同的模型不同的结果
编辑
因此，我能够得出正确的地块，但只有在调整采样间隔后，才能得到以下结果：
t_list = np.linspace(0, 30, 100)
将X打印为：
[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 ...]但这就引出了一个问题:为什么这个系统如此依赖于抽样间隔？
端编辑
我正在尝试重新创建一个简单的matlab系统的微分方程到python使用scipy。我不知道为什么在执行过程中我会收到一个RuntimeWarning: invalid value encountered in double_scalars。我在odeint调用中缺少一个可选参数吗？
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)它印的是：
<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]其中，如下matlab代码所产生的结果是连续的：
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];
endprint语句返回：
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,回答 1
Stack Overflow用户
发布于 2021-09-22 09:21:00
您所处理的ODE系统可能很僵硬。您所遇到的RuntimeWarning是由平方根操作引发的，因为y0的元素在集成过程中变为负值。这是因为积分器的时间步长太大，而且您正在使用的求解器并不适合潜在的刚性系统。增加通过t提供给t_list的元素数量可能会减少时间步长，因此，可以找到解决办法。为了更好地理解正在发生的事情，我鼓励您使用下面的代码片段，它使用了SciPy推荐的更新的SciPy API。特别令人感兴趣的是关键字参数method和max_step。如果使用RK23、RK45、DOP853和LSODA作为一种方法，就会导致解估计结果不佳，因为所有这些求解者都是从非刚性方法开始的。LSODA应该检测刚度并切换到刚性积分器，但随着时间的快速增加，它无法做到这一点。对于所有这些方法，设置max_step=0.5允许以潜在的计算时间代价来处理ODE系统。或者，使用Radau或BDF都会开箱即用，因为这些求解者可以处理僵硬的ODE系统。然而，建议手动提供系统的雅可比，否则它是有限差分近似。
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)https://stackoverflow.com/questions/69244283
复制相似问题
腾讯云开发者
