Немного о картинках

Для поста в проекте N+1 про неустойчивость многозвенных систем управления я решил нарисовать пару картинок. Первая из них была совсем простой: иллюстрация к фразе «любая «хорошая», гладкая функция похожа на линейную в окрестности фиксированной точки». Нужно было взять какую-нибудь функцию и нарисовать её график несколько раз в разных масштабах. В качестве подопытной функции путём проб и ошибок был выбран многочлен $-x-0.45x^2+0.15x^3$.

In [118]:
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
In [119]:
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!

In [120]:
%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$ — собственные значения матрицы, получающейся в модели трёхзвенной системы управления.

In [121]:
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$ (задающиеся экспонентами).

In [122]:
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")

Потом я вошёл во вкус и решил нарисовать фазовый трёхмерной системы. Проблема в том, что рисовать картинки в трёхмерном фазовом пространстве на двумерных (всё ещё) экранах не очень удобно — на них обычно ничего невозможно разглядеть. Но если сделать анимацию, то можно.

In [124]:
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
Out[124]:


Once Loop Reflect