Как построить график, полученный после использования решения resolviv из пакета scipy, для решения системы дифференциальных уравнений в python

Question

Я использую solve_ivp в python, чтобы решить набор дифференциальных уравнений в форме пространства состояний. Мой код выглядит следующим образом:

import numpy as np
import matplotlib.pyplot as plt
from numpy import zeros
from scipy.integrate import solve_ivp
vin=12;
vdon=.7;
X=np.array([[0],[0],[0],[0]])
U= np.array([[vin],[vdon]]);
A1=np.array([[(-.01-10)/1e-5,-1/1e-5,10/1e-5,0] ,[1/1e-8,-1/(.05*1e-8),0,0],[10/1e-8,0,(-10-1)/1e-8,-1/1e-8],[0,0,1/1e-3,0]]);
B1=np.array([[1/1e-5,0],[0,1/(.05*1e-8)],[0,0],[0,0]]);
def conv(t,X):
    xdot= A1.dot(X).reshape(4,1)+B1.dot(U)
    return np.squeeze(np.asarray(xdot))
tspan=[0,.001]

X0=np.array([0,0,0,0])
sol=solve_ivp(conv,tspan,X0)
t=np.linspace(0,.001,10000)
tend=.001
plt.plot(t,X)

Показывает ошибку при использовании нормального plt.plot(t,X) команда. Как построить график между X и t? Пожалуйста, помогите

Answer 1

Немного неясно, что здесь рассчитывается, но, проверяя документацию solve_ivp , кажется, что возвращаемое значение (sol в данном случае) имеет ряд полей. sol.t - список моментов времени. sol.y - это расчетные значения для каждого из 4 компонентов X, соответствующие каждому значению sol.t.

Вот способ построить каждый из 4 компонентов (все они начинаются с 0 потому что X0=(0,0,0,0)). Интерпретация зависит от конкретной решаемой проблемы.

import numpy as np
import matplotlib.pyplot as plt
from numpy import zeros
from scipy.integrate import solve_ivp

vin = 12
vdon = .7
X = np.array([[0], [0], [0], [0]])
U = np.array([[vin], [vdon]])
A1 = np.array([[(-.01 - 10) / 1e-5, -1 / 1e-5, 10 / 1e-5, 0],
               [1 / 1e-8, -1 / (.05 * 1e-8), 0, 0],
               [10 / 1e-8, 0, (-10 - 1) / 1e-8, -1 / 1e-8],
               [0, 0, 1 / 1e-3, 0]])
B1 = np.array([[1 / 1e-5, 0],
               [0, 1 / (.05 * 1e-8)],
               [0, 0],
               [0, 0]])

def conv(t, X):
    xdot = A1.dot(X).reshape(4, 1) + B1.dot(U)
    return np.squeeze(np.asarray(xdot))

tspan = [0, .001]

X0 = np.array([0, 0, 0, 0])
sol = solve_ivp(conv, tspan, X0)
for i in range(sol.y.shape[0]):
    plt.plot(sol.t, sol.y[i], label=f'$X_{i}(t)$')
plt.xlabel('$t$') # the horizontal axis represents the time 
plt.legend() # show how the colors correspond to the components of X
plt.show()

Answer 2

Я нашел более простой способ наметить решение. time_array=sol.t solution_array=sol.y Тогда просто сюжет plt.plot(time_array, solution_array[0,:])