Немного о картинках
Для поста в проекте N+1 про неустойчивость многозвенных систем управления я решил нарисовать пару картинок. Первая из них была совсем простой: иллюстрация к фразе «любая «хорошая», гладкая функция похожа на линейную в окрестности фиксированной точки». Нужно было взять какую-нибудь функцию и нарисовать её график несколько раз в разных масштабах. В качестве подопытной функции путём проб и ошибок был выбран многочлен $-x-0.45x^2+0.15x^3$.
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
oldlims=None
allims=[(-3,3),(-1,1), (-0.3,0.3),(-0.1,0.1)] # область определения для каждого графика
plt.figure(figsize=(11,11./4))
for lims,nextlims,sp in zip(allims,allims[1:]+[None],[141,142,143,144]):
plt.subplot(sp) # нам нужно несколько графиков рядом
X=np.linspace(*lims)
plt.plot(lims,[0,0],color='black') # координатные
plt.plot([0,0],lims,color='black') # оси
plt.plot(X,-X-0.45*X**2+0.15*X**3, color='blue')
plt.plot(X,-X, color='red')
plt.ylim(*lims)
plt.xlim(*lims)
if nextlims:
xmin,xmax=nextlims
plt.plot(*zip([xmin,xmin],[xmin,xmax],[xmax,xmax],[xmax,xmin],[xmin,xmin]),color='black',lw=0.8)
# здесь мы рисуем квадратик, показывающий, какую часть графика увеличиваем
oldlims=lims
Картинка выглядит неплохо, но довольно скучно. Чтобы сделать её повеселее, я решил стилизовать её под рисунок «от руки». Сначала я попробовал нарисовать то же самое действительно вручную, с помощью планшета, но быстро выяснилось, что программирую я всё-таки лучше, чем рисую. Ко всему прочему, на задворках памяти у меня было записано, что я где-то видел пример использовать Python для рисования картинок в стиле комиксов XKCD. Быстрый поиск показал, что всё ещё проще, и сейчас стилизация под XKCD встроена в Matplotlib!
%config InlineBackend.figure_format = 'retina'
with plt.xkcd():
oldlims=None
allims=[(-3,3),(-1,1), (-0.3,0.3),(-0.1,0.1)]
plt.figure(figsize=(11,11./4))
for lims,nextlims,sp in zip(allims,allims[1:]+[None],[141,142,143,144]):
plt.subplot(sp)
X=np.linspace(*lims)
plt.plot(lims,[0,0],color='black',lw=0.7)
plt.plot([0,0],lims,color='black',lw=0.7)
plt.plot(X,-X-0.45*X**2+0.15*X**3)
plt.plot(X,-X)
plt.ylim(*lims)
plt.xlim(*lims)
if nextlims:
xmin,xmax=nextlims
plt.plot(*zip([xmin,xmin],[xmin,xmax],[xmax,xmax],[xmax,xmin],[xmin,xmin]),color='black',lw=0.8)
oldlims=lims
plt.tight_layout(w_pad=0)
plt.savefig("zoom.png")
Так гораздо симпатичнее!
Справившись с этой задачей, я нарисовал ещё комплексные корни третьей степени из числа $-1$ — собственные значения матрицы, получающейся в модели трёхзвенной системы управления.
with plt.xkcd():
X=np.linspace(-2,2)
fig=plt.figure(figsize=(5,5))
ax=fig.add_subplot(111)
ax.plot([-1.5,1.5],[0,0],color='black',lw=0.7)
ax.plot([0,0],[-1.5,1.5],color='black',lw=0.7)
Phi=np.linspace(0,2*np.pi)
ax.plot(np.sin(Phi),np.cos(Phi), lw=0.7)
ax.plot([-1],[0],marker='o',color='green',ms=8,label=u'Устойчивое направление')
ax.plot(*zip([0.5,np.sqrt(3)/2],[0.5,-np.sqrt(3)/2]),marker='o',color='red',linestyle='none',ms=8)
plt.xlim(-1.5,1.5)
plt.ylim(-1.5,1.5)
ax.text(1.45,0.05,u"Вещественная часть",horizontalalignment="right", fontsize=11)
ax.text(-0.05,1.45,u"Мнимая часть",horizontalalignment="right", verticalalignment="top", fontsize=11)
plt.savefig("eigenvalues.png")
А вот так выглядят интегральные кривые (графики решений) уравнения $\dot x=-x$ (задающиеся экспонентами).
with plt.xkcd():
fig=plt.figure(figsize=(12,4))
X=np.linspace(0,3.5)
x0=0
for x0 in range(-2,3,1):
plt.plot(X,x0*np.exp(-X), color='blue' if x0!=0 else 'red')
plt.xlabel(u"время",fontsize=13)
plt.savefig("onedim.png")
Потом я вошёл во вкус и решил нарисовать фазовый трёхмерной системы. Проблема в том, что рисовать картинки в трёхмерном фазовом пространстве на двумерных (всё ещё) экранах не очень удобно — на них обычно ничего невозможно разглядеть. Но если сделать анимацию, то можно.
from matplotlib import animation
from JSAnimation import IPython_display
# Чтобы импортировать JSAnimation, надо его сперва скачать с https://github.com/jakevdp/JSAnimation
# и установить, выполнив команды python setup.py build && python setup.py install
from scipy.integrate import odeint
# эта функция интегрирует дифференциальные уравнения
def F(X,t):
return np.array([X[1],X[2],-X[0]])
# это наша система
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111, projection='3d',xlim=(-3,3),ylim=(-3,3),zlim=(-3,3))
t0=-1.5
t1=4.
for i in xrange(8):
init=np.array([np.sin(i*2*np.pi/8),np.cos(i*2*np.pi/8),0])
ax.plot(*zip(*odeint(F,init,np.linspace(0,t1))),color='blue')
ax.plot(*zip(*odeint(lambda X,t: -F(X,t),init,np.linspace(0,-t0))),color='blue')
# для больше выразительности мы строим траекторию из выбранной точки вперёд и назад во времени
# чтобы постриоть участок «назад во времени» приходится вместо нашего векторного поля брать
# такое же с обратным знаком
x0=odeint(lambda X,t: -F(X,t),init,[0,1.5])[1]
totalframes=16
trajectory=odeint(F,x0,np.linspace(t0,t1,totalframes))
point, = ax.plot([], [], [], marker='o', color='red')
# нашли точку в конце одной из траекторий
def init():
ax.view_init(elev=45, azim=90)
point.set_data(trajectory[0,0:2])
point.set_3d_properties(trajectory[0,2])
return point,
def animate(i):
point.set_data(trajectory[i,0:2])
point.set_3d_properties(trajectory[i,2])
ax.view_init(elev=45,azim=90+i*(360.-90.)/totalframes)
return point,
ani=animation.FuncAnimation(fig, animate, init_func=init,
frames=totalframes, interval=30*360./totalframes)
ani
На получающемся видео хорошо видно устойчивое собственное направление, вдоль которого траектории приближаются к особой точке, и неустойчивая плоскость, на которой траектории убегают от нуля по спиралям. Остаётся сохранить анимацию в файл.
ani.save("animation.mp4", writer='ffmeg', bitrate=1800)
Вот и всё!
Комментарии