import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # function that returns dz/dt def model(z,t,u): x = z[0] y = z[1] dxdt = (-x + u)/2.0 dydt = (-y + x)/5.0 dzdt = [dxdt,dydt] return dzdt # initial condition z0 = [0,0] # number of time points n = 401 # time points t = np.linspace(0,40,n) # step input u = np.zeros(n) # change to 2.0 at time = 5.0 u[51:] = 2.0 # store solution x = np.empty_like(t) y = np.empty_like(t) # record initial conditions x[0] = z0[0] y[0] = z0[1] # solve ODE for i in range(1,n): # span for next time step tspan = [t[i-1],t[i]] # solve for next step z = odeint(model,z0,tspan,args=(u[i],)) # store solution for plotting x[i] = z[1][0] y[i] = z[1][1] # next initial condition z0 = z[1] # plot results plt.plot(t,u,'g:',label='u(t)') plt.plot(t,x,'b-',label='x(t)') plt.plot(t,y,'r--',label='y(t)') plt.ylabel('values') plt.xlabel('time') plt.legend(loc='best') plt.show()
转载于:https://www.cnblogs.com/conpi/p/11102730.html