푸앵카레 단면도 그리기 방법? (Duffing Oscillator)

Question

Duffing 방정식에 대한 간단한 한계 사이클을 성공적으로 보여주는 프로그램을 작성했습니다. 그러나 이제는이 경우 Poincaré 섹션을 그려야합니다.

t*omega = 2*pi*n과 같은 일정한 시간 간격으로 위상 공간 다이어그램의 스냅 샷을 찍어서이 작업을 수행해야합니다. 이 경우에 오메가가 1로 설정되어 있기 때문에 이것은 단지 t = 2*pi*n 일 때입니다. 나는 이것을 시도했지만 Poincaré 섹션을 기대하지는 않습니다.

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은 시스템 기간 중 하나이고 코드는 T=2*p입니다. 실제로 계산하는 것은 (y(n*h), z(n*h))입니다. 여기서 h=0.1은 RK2의 단일 단계의 단계 크기입니다. 빨간색 사각형이

: 다음 그림과 같은 무언가를 얻어야한다

또한 k2y의 인수는 수정 통합으로 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)