ポアンカレ区をプロットする方法は？ （Duffing Oscillator）

Question

Duffing方程式の単純なリミットサイクルを成功裏に示すプログラムを作成しました。しかし、私は今、このケースのためにポアンカレのセクションをプロットする必要があります。

フェーズスペースダイアグラムのスナップショットを一定の時間間隔（例えば、t*omega = 2*pi*n）で取得することで、これを行う必要があります。この場合、オメガが1に設定されているので、これはちょうどt = 2*pi*nのときです。私はこれを試みましたが、私が期待しているポアンカレセクションを得ていません。

はここに私のコードです：

：

program rungekutta 
implicit none 
integer, parameter :: dp = selected_real_kind(15,300) 
integer :: i, n 
real(kind=dp) z, y, t, A, C, D, B, omega, h 
open(unit=100, file="rungekutta.dat",status='replace') 
n = 0 
!constants 
omega = 1.0_dp 
A = 0.25_dp 
B = 1.0_dp 
C = 0.1_dp 
D = 1.0_dp 
y = 0.0_dp 
z = 0.0_dp 
t = 0.0_dp 
do i=1,1000 
call rk2(z, y, t, n) 
n = n + 1.0_dp 
write(100,*) y, z 
end do 

contains 
subroutine rk2(z, y, t, n) !subroutine to implement runge-kutta algorithm 
implicit none 
integer, parameter :: dp = selected_real_kind(15,300) 
integer, intent(inout) :: n 
real(kind=dp) :: k1y, k1z, k2y, k2z, y, z, t, pi 
pi = 4.0*ATAN(1.0) 
h = 0.1_dp 
t = n*2*pi 
k1y = dydt(y,z,t)*h 
k1z = dzdt(y,z,t)*h 
k2z = dzdt(y + (0.5_dp*k1y), z + (0.5_dp*k1z), t + (0.5_dp*h))*h 
k2y = dydt(y, z +(0.5_dp*k1z), t)*h 
y = y + k2y 
z = z + k2z 
end subroutine 

!2nd order ODE split into 2 for coupled Runge-Kutta, useful to define 2 
functions 
function dzdt(y,z,t) 
real(kind=dp) :: y, z, t, dzdt 
dzdt = -A*y**3.0_dp + B*y - C*z + D*sin(omega*t) 
end function 

function dydt(y,z,t) 
real(kind=dp) :: z, dydt, y, t 
dydt = z 
end function 
end program 

私も自分のポアンカレ断面がどのように見えるかの画像を添付しています。

これはx軸とy軸のyです。

そして、私が期待するものの画像：あなたrk2ルーチンで

Answer 1

あなたはステップ長0.1ののステップを実行します。したがって、プロットは、その解像度での解の完全な軌跡です。しかし、その意図は全期間にわたり統合するように思われる。これはそのルーチンのループを必要とするでしょう。

つまり、(y(n*T), z(n*T))のプロットがあります。ここで、Tは、システムの期間の1つです。コードはT=2*pです。実際に計算するのは(y(n*h), z(n*h))です。h=0.1はRK2の1ステップのステップサイズです。赤い四角がある

：あなたは、以下の画像のようなものを取得する必要

もk2yの引数を修正Integratorでuser5713492

のコメントに従って修正する必要がありますポイントはt=n*2*piです。解曲線上の点による指示されたステップサイズは同じh=0.1であり、積分はt=0..300以上である。

def RK2(f,u,times,subdiv = 1): 
    uout = np.zeros((len(times),)+u.shape) 
    uout[0] = u; 
    for k in range(len(times)-1): 
     t = times[k] 
     h = (times[k+1]-times[k])/subdiv 
     for j in range(subdiv): 
      k1 = f(u,t)*h 
      k2 = f(u+0.5*k1, t+0.5*h)*h 
      u, t = u+k2, t+h 
     uout[k+1]=u 
    return uout 

def plotphase(A,B,C,D): 
    def derivs(u,t): y,z = u; return np.array([ z, -A*y**3 + B*y - C*z + D*np.sin(t) ]) 
    N=60 
    u0 = np.array([0.0, 0.0]) 
    t = np.arange(0,300,2*np.pi/N); 
    u = RK2(derivs, u0, t, subdiv = 10) 
    plt.plot(u[:-2*N,0],u[:-2*N,1],'.--y', u[-2*N:,0],u[-2*N:,1], '.-b', lw=0.5, ms=2); 
    plt.plot(u[::N,0],u[::N,1],'rs', ms=4); plt.grid(); plt.show() 

plotphase(0.25, 1.0, 0.1, 1.0)