2016-12-12 45 views
1
from __future__ import division 
import numpy as np 
import matplotlib.pyplot as plt 
import math 
from scipy.integrate import odeint 
from scipy.fftpack import fft, ifft 

def pend(y, t, a, b, ohm): 
    theta, omega, phi = y 
    dydt = [omega, -b*omega-np.sin(theta)-a*np.cos(phi), ohm] 
    return dydt 

b = 1.0/2.0  #beta 
ohm = 2.0/3.0  #capital Omega 
period = 2.0*math.pi/ohm  #driving period 

t0 = 0.0  #initial time 
t = np.linspace(t0,t0+period*10**3,10**3+1)  #time for Poincare map 

theta0 = 0.75 
omega0 = 1.6 
phi0 = 0.8 
y0 = [theta0,omega0,phi0]  #initial conditions 
N = 100       #number of transient points to delete 

a_array = np.linspace(0,1.15,50) #varying parameter of a values 

for a in a_array: 
    sol = odeint(pend,y0,t,args=(a,b,ohm))  #numerical integration of differential equation 
    sol = sol[N:10**3-N]  #removing transients 
    w = sol[:,1]    #frequency 
    A = np.full(len(w),a)  #array of a-values 
    plt.plot(A, w) 
    plt.draw() 

我試圖構造一個當前的分岔圖。在我們使用的方程組中,a是控制參數,我們正在繪製x軸上0到1.15之間的值與特定值a的值(稱爲w)數組之間的關係。我不確定如何從這樣的for循環中繪製事物。我聽說subplots是最好的方法,但我不熟悉實施,可以使用一些幫助。謝謝!在同一圖上繪製多個子圖

+0

現在正在運行的代碼。花點時間跑步。同時,我通常在循環外部移動plt.draw()或plt.show()。 –

回答

0

取消壓縮最後一條命令對我有用。

from __future__ import division 
import numpy as np 
import matplotlib.pyplot as plt 
import math 
from scipy.integrate import odeint 
from scipy.fftpack import fft, ifft 

def pend(y, t, a, b, ohm): 
    theta, omega, phi = y 
    dydt = [omega, -b*omega-np.sin(theta)-a*np.cos(phi), ohm] 
    return dydt 

b = 1.0/2.0  #beta 
ohm = 2.0/3.0  #capital Omega 
period = 2.0*math.pi/ohm  #driving period 

t0 = 0.0  #initial time 
t = np.linspace(t0,t0+period*10**3,10**3+1)  #time for Poincare map 

theta0 = 0.75 
omega0 = 1.6 
phi0 = 0.8 
y0 = [theta0,omega0,phi0]  #initial conditions 
N = 100       #number of transient points to delete 

a_array = np.linspace(0,1.15,50) #varying parameter of a values 

for a in a_array: 
    sol = odeint(pend,y0,t,args=(a,b,ohm))  #numerical integration of differential equation 
    sol = sol[N:10**3-N]  #removing transients 
    w = sol[:,1]    #frequency 
    A = np.full(len(w),a)  #array of a-values 
    plt.plot(A, w) 
plt.show() 

enter image description here

+0

感謝您的及時迴應。我也來到這個情節,但幾乎立即拋棄它,因爲它似乎只顯示1.05到1.15之間的值,這看起來不正確。 – infinitylord

+0

這是你在找什麼? –

+0

只要這是爲a_array中的某個特定值生成w數組,然後將結果繪製在同一個圖上以獲得每個a的值,那麼這就是我正在尋找的 – infinitylord

相關問題