Python으로 Phase Portrait 그리기

6.1.8.py

Strogatz의 Nonlinear Dynamics 교과서에서 6장 연습문제에 보면 van der Pol oscillator가 나온다. 식은 다음과 같이 생겼다. 

이 미분방정식들의 Vector Field를 파이썬으로 그릴 수 있는데 최대한 matplotlib함수를 쓰지 않고(mgrid 등) 그리는 코드를 짰다. 거기에 덧붙여 임의의 점 하나를 잡고 trajectory를 시간에 대해 그리도록 하는 기능을 같이 넣었다. portrait 아래에 나오는 두개의 그래프는 x-t 그리고 y-t 그래프이다. 

trajectory를 따라가보면 시간이 조금 지난뒤에는 진동하는 모습을 보여주며 phase portrait에서도 볼 수 있다. 

Vector plot을 하게되면 trajectory가 뚜렷히 보이지 않는 단점이 있다. vector의 크기가 작을 경우 화살표가 상대적으로 작게 나와 잘 보이지 않기 때문인데, mathematica에서는 StreamPlot이라는 명령어로 trajectory들의 portrait를 살펴보기 쉽다. 하지만 matplotlib에서는 Stream을 plot해주는 함수가 없기 때문에 누군가 만들어 놓은 것이 있다. 아래 링크를 살펴보시면 된다.

http://stackoverflow.com/questions/8296617/how-to-plot-a-streamlines-when-i-know-u-and-v-components-of-velocitynumpy-2d

아래 코드에서간단히 1개의 trajectory를 그렸는데 이를 격자형태로 모두 계산해 그리면 Stream plot이 되는건 뻔한 이야기. -.- (하지만 이건 하고싶지 않아져서 패스)

프로그램을 실행시키면 아래와 같다. (이런걸 왜 넣고 있지.)

참! 예쁘게 나올 필요는 없어서 여백조정은 패스~(사실 이걸 알아보려고 또 20분남짓 소비하는게 싫어서...)

코드는 아래에~!

import numpy as np

import matplotlib.pyplot as plt

def xftn(x, y):

return y

def yftn(x, y):

tmp = -x + y*(1-x*x)

return tmp

def xRK4(x, y, dt):

    k1 = dt * xftn(x, y)

    k2 = dt * xftn(x + 0.5*k1, y + 0.5*k1)

    k3 = dt * xftn(x + 0.5*k2, y + 0.5*k2)

    k4 = dt * xftn(x + k3, y + k3)

    result = x + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6

    return result

def yRK4(x, y, dt):

    k1 = dt * yftn(x, y)

    k2 = dt * yftn(x + 0.5*k1, y + 0.5*k1)

    k3 = dt * yftn(x + 0.5*k2, y + 0.5*k2)

    k4 = dt * yftn(x + k3, y + k3)

    result = y + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6

    return result

Step = 0.1

xLeft = -2

xRight = 2

yTop = 2

yBottom = -2

SIZE = (int)((xRight-xLeft)/Step) + 1  #We only think square plane.

print("Step is", Step)

x = np.empty((SIZE,SIZE), float)

y = np.empty((SIZE,SIZE), float)

u = np.empty((SIZE,SIZE), float)

v = np.empty((SIZE,SIZE), float)

for i in range(SIZE):

for j in range(SIZE):

x[i][j] = xLeft + Step*i

y[i][j] = yBottom + Step*j

u[i][j] = xftn(x[i][j], y[i][j])

v[i][j] = yftn(x[i][j], y[i][j])

xPointList = [0.5]

yPointList = [0.5]

timeList   = [0]

dt = 0.01

for i in range(1500):

xPointList.append(xRK4(xPointList[i], yPointList[i], dt))

yPointList.append(yRK4(xPointList[i], yPointList[i], dt))

timeList.append(timeList[i]+dt)

ax1 = plt.subplot2grid((7,3), (0,0), rowspan=5, colspan=5)

ax1.quiver(x,y,u,v, angles='xy', pivot='middle')

ax1.plot(xPointList, yPointList)

plt.xlim(-2, 2)

plt.ylim(-2, 2)

ax2 = plt.subplot2grid((7,3), (5,0), colspan=3)

ax2.plot(timeList, xPointList)

ax3 = plt.subplot2grid((7,3), (6,0), colspan=3)

ax3.plot(timeList, yPointList)

plt.show()