진자의 방정식을 풀기 위해 Python 스크립트 w 2nd Order Runge Kutta 방법을 사용하여 KE, PE 및 TE의 계산을 어떻게 추가합니까?

Question

KE, PE 및 TE 플롯 스크립트를 프로그래밍하는 방법을 알아낼 수 없습니다. 나는 문제가있는 것 같은 느낌을주는 코드의 부분에 여러 개의 ##########을 포함시켰다.

def pendulum_runge_kutta(theta0,omega0,g,tfinal,dt): 
    # initialize arrays 
    t = np.arange(0.,tfinal+dt,dt)    # time array t 
    npoints = len(t) 
    theta = np.zeros(npoints)     # position array theta 
    omega = np.zeros(npoints)     # position array omega 
    Ke = np.zeros(npoints) 
    Pe = np.zeros(npoints) 
    L=4 
    g = 9.81 
    t2=np.linspace(0,tfinal,1000) 
    theta0=0.01 
    omega0=0 

    # exact solution for 
    thetaExact = theta0*np.cos((g/L)**(1/2)*t2) 

    # SECOND ORDER RUNGE_KUTTA SOLUTION 
    theta[0] = theta0 
    omega[0] = omega0 
    #Ke[0] = ####################################### 
    #Pe[0] =###################################### 
    m=1.0 
    for i in range(npoints-1): 
     # compute midpoint position (not used!) and velocity   
     thetamid = theta[i] + omega[i]*dt 
     omegamid = omega[i] - (g/L)*np.sin(theta[i])*dt/2 

     # use midpoint velocity to advance position    
     theta[i+1] = theta[i] + omegamid*dt 
     omega[i+1] = omega[i] -(g/L)*np.sin(thetamid)*dt/2 

     ###########calculate Ke, Pe, Te############ 
     Ke[i+1] = 0.5*m*(omega[i+1]*L)**2 
     Pe[i+1] = m*g*L*np.sin(theta[i+1]) 
    Te = Ke+Pe 

    #plot result of Ke, Pe, Te 
    pl.figure(1) 
    pl.plot(t,Ke,'c-',label='kinetic energy') 
    pl.plot(t,Pe,'m-',label='potential energy') 
    pl.plot(t,Te,'g-',label='total energy') 
    pl.title('Ke, Pe, and Te') 
    pl.legend(loc='lower right') 
    pl.show() 

    #now plot the results 
    pl.figure(2) 
    pl.plot(t,theta,'ro',label='2oRK') 
    pl.plot(t2,thetaExact,'r',label='Exact') 
    pl.grid('on') 
    pl.legend() 
    pl.xlabel('Time (s)') 
    pl.ylabel('Theta') 
    pl.title('Theta vs Time') 
    pl.show() 

#call function  
pendulum_runge_kutta(0.01, 0, 9.8, 12, .1) 

Answer 1

중간 점 방법의 공식은 1 차 시스템의 모든 구성 요소에 균일하게 적용됩니다. 따라서 중간 점 값은 모두 dt/2이고 전체 시간 단계는 전체 dt이어야합니다.

잠재 에너지는 sin(x), C-cos(x)의 정수를 포함해야합니다. 예를 들어, 시점 i+1에서 에너지에 대한

Pe[i+1] = m*g*L*(1-np.cos(theta[i+1])) 

수식은 시점 0에 적용됩니다. 예를 들어 길이가 L=4 인 경우 대량 선언을 위로 이동해야합니다.

최종 설명에서 표시된 정확한 해결책은 작은 각도 근사를위한 것입니다. 일반 육체 진자에 대해 정확하게 표현할 수있는 정확한 해법이 없습니다. 주어진 진폭이 theta0=0.01 인 경우, 작은 각 근사가 충분해야하지만 큰 스윙의 경우이를 알고 있어야합니다.