Блог

Главная - Вопросы по программированию - используя скрипт на Python с методом Рунге-Кутты 2-го порядка для решения уравнения маятника, как мне добавить вычисление KE, PE и TE?

используя скрипт на Python с методом Рунге-Кутты 2-го порядка для решения уравнения маятника, как мне добавить вычисление KE, PE и TE?

#python #physics #runge-kutta

#python #физика #рунге-Кутта

Вопрос:

Я не могу понять, как запрограммировать мой скрипт для построения графиков 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)

Комментарии:

1. Не могли бы вы, пожалуйста, указать ссылку на математические уравнения, которые вы пытаетесь решить? Является ли единственной проблемой присвоение Ke [0] и Pe [0]?

2. @Jalo: Обычный физический маятник, m*L*theta'' m*g*sin(theta)=0 который при малых углах может быть аппроксимирован уравнением гармонических колебаний theta'' g/L*theta=0 .

Ответ №1:

Формулы для метода средней точки применяются равномерно ко всем компонентам системы первого порядка. Таким образом, это должно быть dt/2 как для средних значений, так и для полных dt для полного временного шага.

Потенциальная энергия должна содержать интеграл от sin(x) , C-cos(x) . Например, 

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

Формулы для энергий в момент времени i 1 также применяются в момент времени 0 . Вам нужно переместить объявление массы вверх, например, после строки длины L=4 .

И последнее замечание, указанное точное решение предназначено для приближения с малым углом, для общего физического маятника нет конечно выражаемого точного решения. Для заданной амплитуды theta0=0.01 приближение малого угла должно быть достаточно хорошим, но имейте в виду это для больших колебаний.