Маневр перехода с плоским вращением: почему угловой момент непостоянен?

Краткий ответ: это проблема решателя ODE.

При настройке допусков относительной и абсолютной погрешности достигаются лучшие результаты. Например, optionsпеременная в коде MATLAB может быть определена следующим образом:

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

Кроме того, чтобы начальный угловой момент был 3000вместо 3000.0045(как было указано в этом ответе ), уравнение для расчета углового момента можно решить для начальной скорости вращения малой оси следующим образом:

omega_z = ( ( 3000^2 - ( Ix * omega_x )^2 ) / ( Iz^2 ) )^( 1/2 );

Наконец, чтобы настроить оси графика углового момента, добавьте:

axis([0 1000 2999 3001])

Получается следующая фигура:

Я переписал ваш матлаб на Python (что действительно легко сделать, и это не совпадение!) И использовал интегратор SciPy ODE по умолчанию. Проверяя, info['mused']я вижу, что он всегда использовал алгоритм Адамса (нежесткий), поэтому проблема не в жесткости.

Почему бы просто не сбросить «старый и сломанный» Matlab и не принять «новый» Python? ( ссылка на MIB )

Я предполагаю, что вам следует начать читать документацию по Matlab ode45и снизить допуск к ошибкам. Я предполагаю, что это ode45то, что я бы назвал RK45, и это переменный размер шага 4/5 порядка Рунге-Кутта. Это должно быть хорошо, но если вы не укажете допуск, он выберет его для вас.

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

Короче говоря, у меня он работает просто отлично.

Начальный угловой момент не равен 3000 с использованием значений инициализации, но на самом деле это 3000.0045то место, где начинается симуляция. Он дрейфует 3000.001через 1000 секунд. Это изменение 0,003. На графике показано отклонение от четных 3000 (или 3e3).

Я могу немного улучшить это, чтобы общий дрейф составлял всего около 0,001, вставив rtol=1E-10или rtol=1E-11в вызове ODEint.

def omega_dynamics(x_0, t, I, mu, T, J):

    # % Initial state vector
    omega_x, omega_y, omega_z, sigma_x, sigma_y, sigma_z = x_0[:6]

    # % Constants
    mu_x, mu_y, mu_z = mu[:3]
    Ix, Iy, Iz = I[:3]
    Tx, Ty, Tz = T[:3]

    # % Differential equations
    omega_dot_x = ((Iy - Iz)*omega_y*omega_z + mu_x*sigma_x + Tx)/(Ix - J)
    omega_dot_y = ((Iz - Ix)*omega_z*omega_x + mu_y*sigma_y + Ty)/(Iy - J)
    omega_dot_z = ((Ix - Iy)*omega_x*omega_y + mu_z*sigma_z + Tz)/(Iz - J)
    sigma_dot_x = -omega_dot_x - mu_x*sigma_x/J - omega_y*sigma_z + omega_z*sigma_y
    sigma_dot_y = -omega_dot_y - mu_y*sigma_y/J - omega_z*sigma_x + omega_x*sigma_z
    sigma_dot_z = -omega_dot_z - mu_z*sigma_z/J - omega_x*sigma_y + omega_y*sigma_x

    # % Return vector
    dx = np.array([omega_dot_x, omega_dot_y, omega_dot_z, sigma_dot_x, sigma_dot_y, sigma_dot_z])

    return dx

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %%%%% Example Rahn & Barba (1991) %%%%%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

J        = 18.0;                                             # Principal moments of inertia of spherical propellant slug
I        = np.array([2000,   1500,   1000], dtype=float)   # Principal moments of inertia of spacecraft including slug
omega_0  = np.array([0.1224, 0,      2.99], dtype=float)   # Spacecraft body rates
sigma_0  = np.array([0,      0,      0   ], dtype=float)   # Relative rates between spacecraft body and propellant slug
mu       = np.array([30,     30,     30  ], dtype=float)   # Viscous damping coefficient
T        = np.array([0,      0,      0   ], dtype=float)   # Torques about principal axes

h_init   = np.sqrt(((I * omega_0 + J * sigma_0)**2).sum())  # % Calculate angular momentum

print "h_init: ", h_init

times        = np.linspace(0, 1000, 10001)
X0           = np.hstack((omega_0, sigma_0))

x1, info = ODEint(omega_dynamics, X0, times, args=(I, mu, T, J), full_output=True)
print x1.shape

h          = np.sqrt(((I * x1[:, :3] + J * x1[:,3:])**2).sum(axis=1))  # % Calculate angular momentum

if True:
    plt.figure()

    plt.subplot(2, 2, 1)
    plt.plot(times, x1[:, 2])
    plt.ylim(-2, 3)
    plt.title('Minor Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 2)
    plt.plot(times, x1[:, 1])
    plt.title('Intermediate Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 3)
    plt.plot(times, x1[:, 0])
    plt.title('Major Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 4)
    plt.plot(times, h)
    plt.title('Angular Momentum')
    plt.xlabel('Time - seconds')
    plt.ylabel('Nms')

    plt.show()