Strogatz의 Nonlinear Dynamics 교과서에서 6장 연습문제에 보면 van der Pol oscillator가 나온다. 식은 다음과 같이 생겼다.
이 미분방정식들의 Vector Field를 파이썬으로 그릴 수 있는데 최대한 matplotlib함수를 쓰지 않고(mgrid 등) 그리는 코드를 짰다. 거기에 덧붙여 임의의 점 하나를 잡고 trajectory를 시간에 대해 그리도록 하는 기능을 같이 넣었다. portrait 아래에 나오는 두개의 그래프는 x-t 그리고 y-t 그래프이다.
Vector plot을 하게되면 trajectory가 뚜렷히 보이지 않는 단점이 있다. vector의 크기가 작을 경우 화살표가 상대적으로 작게 나와 잘 보이지 않기 때문인데, mathematica에서는 StreamPlot이라는 명령어로 trajectory들의 portrait를 살펴보기 쉽다. 하지만 matplotlib에서는 Stream을 plot해주는 함수가 없기 때문에 누군가 만들어 놓은 것이 있다. 아래 링크를 살펴보시면 된다.
아래 코드에서간단히 1개의 trajectory를 그렸는데 이를 격자형태로 모두 계산해 그리면 Stream plot이 되는건 뻔한 이야기. -.- (하지만 이건 하고싶지 않아져서 패스)
프로그램을 실행시키면 아래와 같다. (이런걸 왜 넣고 있지.)
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()
6.1.8_variation_uploaded.py
