Using ODE45 to Solve a Differential Equation
Introduction
For this tutorial, I will demonstrate how to use the ordinary differential equation solvers within MATLAB to numerically solve the equations of motion for a satellite orbiting Earth.
For two-body orbital mechanics, the equation of motion for an orbiting object relative to a much heavier central body is modeled as:
Where μ is the gravitational parameter of earth (398600 km3/s2)
The state-space representation of this equation becomes:
![\left[ \begin{matrix} {\dot{\vec{r}}} \\ {\dot{\vec{v}}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{0}_{3\times 3}} & {{\text{I}}_{3\times 3}} \\ {\mu }/{{{r}^{3}}}\; & {{0}_{3\times 3}} \\ \end{matrix} \right]\left[ \begin{matrix} {\vec{r}} \\ {\vec{v}} \\ \end{matrix} \right]](http://s0.wp.com/latex.php?latex=%5Cleft%5B+%5Cbegin%7Bmatrix%7D+%7B%5Cdot%7B%5Cvec%7Br%7D%7D%7D+%5C%5C+%7B%5Cdot%7B%5Cvec%7Bv%7D%7D%7D+%5C%5C+%5Cend%7Bmatrix%7D+%5Cright%5D%3D%5Cleft%5B+%5Cbegin%7Bmatrix%7D+%7B%7B0%7D_%7B3%5Ctimes+3%7D%7D+%26+%7B%7B%5Ctext%7BI%7D%7D_%7B3%5Ctimes+3%7D%7D+%5C%5C+%7B%5Cmu+%7D%2F%7B%7B%7Br%7D%5E%7B3%7D%7D%7D%5C%3B+%26+%7B%7B0%7D_%7B3%5Ctimes+3%7D%7D+%5C%5C+%5Cend%7Bmatrix%7D+%5Cright%5D%5Cleft%5B+%5Cbegin%7Bmatrix%7D+%7B%5Cvec%7Br%7D%7D+%5C%5C+%7B%5Cvec%7Bv%7D%7D+%5C%5C+%5Cend%7Bmatrix%7D+%5Cright%5D&bg=ffffff&fg=000&s=1 "\left[ \begin{matrix} {\dot{\vec{r}}} \\ {\dot{\vec{v}}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{0}_{3\times 3}} & {{\text{I}}_{3\times 3}} \\ {\mu }/{{{r}^{3}}}\; & {{0}_{3\times 3}} \\ \end{matrix} \right]\left[ \begin{matrix} {\vec{r}} \\ {\vec{v}} \\ \end{matrix} \right]")
MATLAB has many ODE solvers available for the coder. The selection of the appropriate solver is dependent on the type of ODE you are solving and the desired accuracy. For this problem, we will use the ode45 solver which uses a Runge-Kutta iterative method to achieve 4th and 5th order accuracy. I recommend that students write their own Runge-Kutta function to better understand this algorithm prior to adopting that MATLAB internal function. The ode45 function within MATLAB uses the Dormand-Prince formulation.
To understand the input parameters for the ode45 function, type “doc ode45” and “doc odeset” in the MATLAB command window.
Now Let’s Get Started
For this problem, the equation of motion for the satellite will be coded as an anonymous function. Anonymous functions are convenient since they can be written within another script of function, but they can only contain one executable statement. For this problem, anonymous functions are an excellent choice to code the state-space equation. NOTE: Any variable referenced by an anonymous function that is not an input to the function, must be define before the function within the code. For example: mu.
% Define Constants
mu = 3.986012E5; % km^3/s^2
% Define the equation of motion
EoM = @(t, x) [zeros(3,3), eye(3); -(mu/norm(x(1:3,1))^3)*eye(3), zeros(3,3)]*x;
% Determine Orbital Period, this is just to determine how long to simulate for
E = norm(V0)^2/2 - mu/norm(R0);
a = -mu/(2*E);
n = sqrt(mu/a^3);
T = 2*pi/n;
% Now lets plot the results
options = odeset('MaxStep', 5); % Set maximum step size as 5 seconds
[time, State] = ode45(EoM, [0 T], X0, options);
How Did I Know How Long to Integrate For?
Simply knowing R and V at any point in time and the gravitational parameter of the central body, the period of the orbit can be determined.
First the total energy of the orbit needs to be calculated. The total energy is composed of the kinetic and potential energies and remains constant throughout the orbit:
The semi-major axis of the orbit can then be determined:
Finally, from the semi-major axis, the orbital period becomes:
NOTE: The orbital period can only be determined for circular and elliptical orbits (E < 0).