概要:
力学系を勉強するときに、平面上にベクトル場(vector fields)を書いて、微分方程式の解を求めるときがある。このときにPythonのmatplotlibを使って、ベクトル場を書く。
のベクトル場
Input:
# Vector fields of the differential equation x'' = -x import matplotlib.pyplot as plt import numpy as np plt.figure() (LX, LY) = (5, 5) gridwidth=0.4 X, Y= np.meshgrid(np.arange(-LX, LX, gridwidth), np.arange(-LY, LY,gridwidth)) R = np.sqrt(X**2+Y**2) X1,Y1=0,0 R1=np.sqrt(X**2+Y**2) plt.plot(X1,Y1,'o',color='blue') #Vector fields F(X, Y) = (U(X, Y), V(X, Y) ) = (Y/R1, -X/R1) U = Y/R1 V = -X/R1 plt.quiver(X,Y,U,V,color='black',angles='xy',scale_units='xy', scale=4.5) plt.xlim([-LX,LX]) plt.ylim([-LY,LY]) # graph plt.grid() plt.draw() plt.show()
Output:
の解曲線
Input:
# Stream plots of the differential equation x'' = -x import numpy as np import matplotlib.pyplot as plt w = 5 Y, X = np.mgrid[-w:w:100j, -w:w:100j] #Vector fields F(X, Y) = (U(X, Y), V(X, Y) ) = (Y, -X) U = Y V = -X #the fixed point at origin seed_points = np.array([[0], [0]]) plt.plot(seed_points[0], seed_points[1], 'bo') #graph plt.streamplot(X, Y, U, V) plt.show()
Output:
微分方程式 について
と定義すると、次のように書き換えることができる。
この方程式で定められるベクトル場を相空間(Phase space)に書けばよい。
のベクトル場
Input:
import matplotlib.pyplot as plt import numpy as np plt.figure() (LX, LY) = (5, 5) gridwidth=0.4 X, Y= np.meshgrid(np.arange(-LX, LX, gridwidth), np.arange(-LY, LY,gridwidth)) R = np.sqrt(X**2+Y**2) #fixed point at origin X1,Y1=0,0 R1=np.sqrt(X**2+Y**2) plt.plot(X1,Y1,'o',color='blue') #Vector fields U = Y V = -np.sin(X) plt.quiver(X,Y,U,V,color='black',angles='xy',scale_units='xy', scale=4.5) plt.xlim([-LX,LX]) plt.ylim([-LY,LY]) # graph plt.grid() plt.draw() plt.show()
Output:
各点のベクトルを正規化されていない。つまり、各点のベクトルの長さが1となっていない。しかしこちらのほうが幾分わかりやすかった。
の解曲線
Input:
import numpy as np import matplotlib.pyplot as plt w = 4 Y, X = np.mgrid[-w:w:1000j, -w:w:1000j] U = Y V = -np.sin(X) plt.streamplot(X, Y, U, V, color='black') seed_points = np.array([[0], [0], [3.05], [-3.15]]) plt.plot(seed_points[0], seed_points[1], 'bo') plt.plot(seed_points[2], seed_points[0], 'bo') plt.plot(seed_points[3], seed_points[0], 'bo') plt.show()
Output:
