Расчет положения спутника с использованием Runge-Kutta 4

Question

Расчет положения спутника с использованием Runge-Kutta 4

Моя проблема связана с методом Рунге-Кутты 4 (RK4) и правильными шагами итерации, необходимыми для вектора состояния орбитального спутника. Приведенный ниже код (на Python) описывает движение на основе описания по этой ссылке ( http://www.navipedia.net/index.php/GLONASS_Satellite_Coordinates_Computation):

    if total_step_number != 0:   
        for i in range(1, total_step_number+1):                             
            #Calculate k1                
            k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
            k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
            k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[2]

            #Intermediate step to bridge k1 to k2
            XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
            XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
            XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
            XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
            XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
            XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
            radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

             ....

Существует больше кода, однако я хочу ограничить то, что я пока показываю, так как это промежуточные шаги, которые я больше всего заинтересован в решении. В основном, для тех, кто знаком с векторами состояний и использует RK4, вы можете видеть, что положение и скорость обновляются на промежуточном этапе, но не ускорение. Мой вопрос связан с расчетом, необходимым для того, чтобы обновить тоже ускорение. Это началось бы:

XYZDDot[0] = ...
XYZDDot[1] = ...
XYZDDot[2] = ...

... но что будет после, не очень понятно. Любые советы приветствуются.

Ниже приведен полный код:

        for j in h_step_values:
            h = j    
            if h > 0:
                one_way_iteration_steps = one_way_iteration_steps -1
            elif h < 0:
                one_way_iteration_steps = one_way_iteration_steps +1
                XYZ = initial_XYZ
            #if total_step_number != 0:   
            for i in range(0, one_way_iteration_steps):                             
                #Calculate k1                
                k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
                k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
                k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k1 to k2
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
                XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
                XYZDDot2[0] = XYZDDot[0] + (k1[0] * h / 2)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
                XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
                XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k2  
                k2[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
                k2[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
                k2[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k2 to k3
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k2[0] * h**2 / 8)
                XYZDot2[0] = XYZDot[0] + (k2[0] * h / 2)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k2[1] * h**2 / 8)
                XYZDot2[1] = XYZDot[1] + (k2[1] * h / 2)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k2[2] * h**2 / 8)
                XYZDot2[2] = XYZDot[2] + (k2[2] * h / 2)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k3  
                k3[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1]) 
                k3[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
                k3[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k3 to k4
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h) + (k3[0] * h**2 / 2)
                XYZDot2[0] = XYZDot[0] + (k3[0] * h)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h) + (k3[1] * h**2 / 2)
                XYZDot2[1] = XYZDot[1] + (k3[1] * h)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h) + (k3[2] * h**2 / 2)
                XYZDot2[2] = XYZDot[2] + (k3[2] * h)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k4 
                k4[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
                k4[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0]) 
                k4[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]


                for p in range(3):
                    XYZ[p] = XYZ[p] + XYZDot[p] * h + h**2 * ((k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 12)
                    XYZDot[p] = XYZDot[p] + (h * (k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 6)

                radius = np.sqrt((XYZ[0])**2 + (XYZ[0])**2 + (XYZ[0])**2)

3

python runge-kutta orbital-mechanics satellite-navigation

Источник

user5257286 19 авг '16 в 09:32

1 ответ

Другие вопросы по тегам python runge-kutta orbital-mechanics satellite-navigation

user3088138 21 авг '16 в 17:15 2016-08-21 17:15 · Answer 1 · 2016-08-21 17:15

Уравнение, которое вы решаете, имеет вид

ddot x = a(x)

где a(x) это ускорение, которое вычисляется в вашем k1 вычисление. Действительно, система первого порядка будет

dot v = a(x)
dot x = v

Таким образом, реализация RK4 начинается с

k1 = a(x)
l1 = v

k2 = a(x+l1*h/2) = a(x+v*h/2)
l2 = v+k1*h/2

и т. д. Использование l1,l2,... кажется неявным в коде, вставляя эти линейные комбинации непосредственно там, где они встречаются.

Короче говоря, вы не пропустите ускорение вычислений, это основная часть фрагмента кода.

Обновление: (8/22) Чтобы приблизиться к цели промежуточных шагов моста, абстрактный код должен быть прочитан (с (* .. *) обозначает комментарии или ненужные вычисления)

k1 = a(x)                    (* l1 = v *)

x2 = x + v*h/2               (* v2 = v + k1*h/2 *)

k2 = a(x2)                   (* l2 = v2 *)

x3 (* = x + l2*h/2 *) 
   = x + v*h/2 + k1*h^2/4    (* v3 = v + k2*h/2 *)

k3 = a(x3)                   (* l3 = v3 *)

x4 (* = x + l3*h *)
   = x + v*h + k2*h^2/2      (* v4 = v + k3*h *)

k4 = a(x4)                   (* l4 = v4 *)


delta_v = ( k1+2*(k2+k3)+k4 ) * h/6

delta_x (* = ( l1+2*(l2+l3)+l4 ) * h/6 *)
        = v*h + (k1+k2+k3) * h^2/6