How to plot the graph obtained after using solve_ivp from scipy package to solve a set of differential equations in python

I’m using solve_ivp in python to solve a set of differential equations in state-space form.My code is as follows:

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)

It shows error while using the normal plt.plot(t,X) command. How do I plot the graph between X and t? Please help

Answer

It’s a bit unclear what’s being calculated here, but checking solve_ivp‘s documentation, it seems the return value (sol in this case) has a number of fields. sol.t is the list of time points. sol.y are the calculated values for each of the 4 components of X, corresponding to each value of sol.t.

Here is a way to plot each of the 4 components (they all start at 0 because X0=(0,0,0,0)). The interpretation depends on the exact problem being solved.

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()