하는 일

대학의 후배가 선형계의 축퇴의 이미지가 솟지 않는다고 했기 때문에 2차원 선형계의 고유 벡터가 융합해 가는 애니메이션을 만들려고 했습니다. 풀고있는 방정식은 매우 간단하며 다음과 같은 형태를 취합니다.

\frac{d}{dt}\Bigl(\begin{matrix}
x\\
y
\end{matrix} \Bigl)=
\Bigl(\begin{matrix}
-1 & -1\\
0 & -1-\epsilon
\end{matrix}\Bigl)\Bigl(\begin{matrix}x\\
y\end{matrix}\Bigl)

$2\times 2$ 행렬의 $(2,2)$ 성분에만 섭동 $\epsilon$를 넣어 두고 축퇴를 미리 풀고 있습니다. 그런 다음 for 문을 돌려 $\epsilon $을 점진적으로 변경하고 각 단계마다 여러 초기 조건에서 scipy의 odeint로 적분합니다.
아래의 샘플에서는 안정적인 노드를 축퇴시키고 있습니다. 좀 더 for문을 돌리면 원점이 안장으로 변화하는 것도 볼 수 있다고 생각하고, 섭동의 넣는 방법에 따라서는 노드와 나선의 전이도 볼 수 있을 것입니다.

실제로 작성한 코드

반지성적인 코드이므로 공개할지 반나절 고민했습니다만, 누군가에게 의미있는 코멘트를 받는 것을 기대해 공개합니다. for 문에서 단계별 섭동 파라미터 $ e (=\epsilon) $를 단계 수에 따라 계산하고, 그것을 바탕으로 미분 방정식을 풀고 있습니다.

Linear_degenerate.ipynb

%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi
import matplotlib.animation as animation
# degeneration of linear system
def func(v, t, e):
    return [ -v[0] - v[1], -(1+e)*v[1] ] #  2-D linear system

#def func(v,t):
#    return [ -v[1] -0.5*v[0]*(v[0]**2+v[1]**2), v[0] -0.5*v[1]*(v[0]**2 + v[1]**2)]

# time length and unit step    
dt = 0.1
T = np.arange(0.0, 20, dt)
v0 = [[17,7], [17,17], [7,17], [-17,-7], [-17,-17], [-7,-17],[-7,17], [7,-17]]

# make animation
fig = plt.figure()

ims = []

for k in range(200):
    e = (150-k)/50
    #each initial state
    space = []
    for v, color in zip(v0, 'rbgrbgkk'):
        u = odeint(func, v, T, args=(e,))
        space += plt.plot(u[:,0], u[:,1], color)

    ims.append(space)


    '''u0 = odeint(func, v0[0], T, args=(e,))
    space0 = plt.plot(u0[:,0], u0[:,1], 'r')
    u1 = odeint(func, v0[1], T, args=(e,))
    space1 = plt.plot(u1[:,0], u1[:,1], 'b')
    u2 = odeint(func, v0[2], T, args=(e,))
    space2 = plt.plot(u2[:,0], u2[:,1], 'g')
    u3 = odeint(func, v0[3], T, args=(e,))
    space3 = plt.plot(u3[:,0], u3[:,1], 'r')
    u4 = odeint(func, v0[4], T, args=(e,))
    space4 = plt.plot(u4[:,0], u4[:,1], 'b')
    u5 = odeint(func, v0[5], T, args=(e,))
    space5 = plt.plot(u5[:,0], u5[:,1], 'g')
    u6 = odeint(func, v0[6], T, args=(e,))
    space6 = plt.plot(u6[:,0], u6[:,1], 'k')
    u7 = odeint(func, v0[7], T, args=(e,))
    space7 = plt.plot(u7[:,0], u7[:,1], 'k')

    ims.append(space0+space1+space2+space3+space4+space5+space6+space7)
    '''
plt.xlim(-17 ,17) 
plt.ylim(-17 ,17)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
ani = animation.ArtistAnimation(fig, ims, interval=27, blit=True, repeat=True)

plt.show()

실행 예

$x$와는 1차 독립이었던 고유 방향이 점차 $x$에 붕괴해 가는 것이 보인다고 생각합니다.

문제점

앞에서 설명한 코드는 문제로 가득합니다.
1 : for 문에서 비슷한 작업을 여러 번 수행합니다.
2:자동화하고 싶었지만 초기값에 따른 색의 변경을 어떻게 할까 생각하고 있지 않았다
3: 모처럼 축퇴를 보고 있는데 고유 방향의 플롯을 잊어버린다
등을 들 수 있습니다.

더 현명한 코드를 작성할 수 있으면 나중에 어떻게든 할 것입니다. 가르쳐 준 부드러운 분 덕분에 문제 1, 2가 해결되었습니다. 다시 한번 감사드립니다. 해결해 주신 문제의 부분입니다만, 일단 자신이 쓴 더러운 코드도 코멘트 아웃으로 남겨두기로 했습니다.
초기 조건의 장소도 좀 더 궁리할 수 있으면, 라고 생각했습니다만 계에 의해 보기 쉬운 장소는 다르기 때문에 결국 손으로 넣을 수밖에 없는 것일까라고.