Блог

Редактировать

Таким образом, я смог создать правильный график, но это было только после настройки интервала выборки на следующий:

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?

питон:

 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 amp;amp;  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
 

инструкция print возвращает:

 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,
 

Система ODE, с которой вы имеете дело, скорее всего, жесткая. Проблема RuntimeWarning , с которой вы сталкиваетесь, возникает в результате операций с квадратным корнем, поскольку элементы y0 становятся отрицательными во время интеграции. Это происходит из-за того, что время работы интегратора слишком велико, а используемый вами решатель плохо подходит для потенциально жестких систем. Увеличение количества элементов, предоставляемых для t сквозного t_list , вероятно, уменьшает время и, следовательно, позволяет обойти проблему. Чтобы лучше понять, что происходит, я рекомендую вам поиграть с приведенным ниже фрагментом, в котором используется более новый solve_ivp API, рекомендованный SciPy. Особый интерес представляют ключевые слова аргументы 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)