using a Python script w 2nd Order Runge Kutta method to solve the equation of a pendulum, how do I add a calculation of the KE, PE, and TE?
<p>I am unable to figure out how to program my script to plot KE, PE, and TE. I have included multiple ########### in the parts of my code where I feel the problem lies.  </p>

<pre><code>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)
</code></pre>

<p>You are missing a factor 1/2 in both theta updates. The formulas for the midpoint method apply uniformly to all components of the first order system.</p>

<p>The potential energy should contain the integral of <code>sin(x)</code>, <code>C-cos(x)</code>. For instance, </p>

<pre><code>Pe[i+1] = m*g*L*(1-np.cos(theta[i+1]))
</code></pre>

<p>The formulas for the energies at time point <code>i+1</code> also apply at time point <code>0</code>.</p>

<p>And finally, the indicated exact solution is for the small angle approximation, there is no finitely expressible exact solution for the general physical pendulum. For the given amplitude of <code>theta0=0.01</code>, the small angle approximation should be good enough, but be aware of this for larger swings.</p>

python physics runge-kutta