Solve Celestial Mechanics Problem with High-Order Solvers

This example shows how to use ode78 and ode89 to solve a celestial mechanics problem that requires high accuracy in each step from the ODE solver for a successful integration. Both ode45 and ode113 are unable to solve the problem using the default error tolerances. Even with more stringent error thresholds, ode89 performs best on the problem due to the high accuracy of the Runge-Kutta formulas it uses in each step.

Problem Description

The Pleiades problem is a celestial mechanics problem describing the gravitational interactions of seven stars [1]. This cluster of stars is also referred to as The Seven Sisters, and it is visible to the human eye in the night sky due to its proximity to Earth [2]. The system of equations describing the motion of the stars in the cluster consists of 14 nonstiff second-order differential equations, which produce a system of 28 equations when rewritten in first-order form.

Newton's second law of motion relates the force applied to a body to its rate of change in momentum over time,

${\mathit{F}}_{\mathit{i}}=\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{p}}_{\mathit{i}}$.

The momentum (${\mathit{p}}_{\mathit{i}}={\mathit{m}}_{\mathit{i}}{\mathit{v}}_{\mathit{i}}$) has separate x- and y-components. At the same time, the universal law of gravitation describes the force working on body i from body j as

${\mathit{F}}_{\mathrm{ij}}=\mathit{g}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\Vert {\mathit{p}}_{\mathit{i}}-{\mathit{p}}_{\mathit{j}}\Vert }^2}{\mathit{d}}_{\mathrm{ij}}$.

The term ${\mathit{d}}_{\mathrm{ij}}=\frac{{\mathit{p}}_{\mathit{j}}-{\mathit{p}}_{\mathit{i}}}{\Vert {\mathit{p}}_{\mathit{j}}-{\mathit{p}}_{\mathit{i}}\Vert }$ gives the direction of the distance between the bodies, and the masses of the bodies are ${\mathit{m}}_{\mathit{i}}=\mathit{i}$ for $\mathit{i}=1,2,...,7$. For a system with many bodies, the force applied to any individual body is the sum of its interactions with all others, so

$${\mathit{F}}_{\mathit{i}}=\sum _{\mathit{i}\ne \mathit{j}}{\mathit{F}}_{\mathrm{ij}}.$$

Setting the gravitational constant $\mathit{g}$ equal to 1 and solving yields a system of second-order equations describing the evolution of the x- and y-components over time.

$${{{\mathit{x}}_{\mathit{i}}}^{\prime }}^{\prime }={\mathit{f}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)=\sum _{\mathit{j}\ne \mathit{i}}\frac{{\mathit{m}}_{\mathit{j}}\left({\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}}\right)}{{\mathit{r}}_{\mathrm{ij}}^{3/2}},$$

$${{{\mathit{y}}_{\mathit{i}}}^{\prime }}^{\prime }={\mathit{h}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)=\sum _{\mathit{j}\ne \mathit{i}}\frac{{\mathit{m}}_{\mathit{j}}\left({\mathit{y}}_{\mathit{j}}-{\mathit{y}}_{\mathit{i}}\right)}{{\mathit{r}}_{\mathrm{ij}}^{3/2}},$$

where ${\mathit{r}}_{\mathrm{ij}}=\sqrt{{\left({\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{j}}\right)}^2+{\left({\mathit{y}}_{\mathit{i}}-{\mathit{y}}_{\mathit{j}}\right)}^2}$. Because these two equations apply for each of the seven stars in the system, 14 second-order differential equations ($\mathit{i}=1,2,...,7)$ describe the entire system.

The initial conditions for the system, as given in [1], are:

$$\{\begin{array}{ll}{x}_0& =(3,3,-1,-3,2,-2,2)\\ y_0& =(3,-3,2,0,0,-4,4)\\ x_0^{\prime }& =(0,0,0,0,0,1.75,-1.5)\\ y_0^{\prime }& =(0,0,0,-1.25,1,0,0)\end{array}$$

To solve this system of ODEs in MATLAB®, you must code the equations into a function before calling the solvers ode78 and ode89. You can either include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Code Equations

The ODE solvers in MATLAB require equations to be written in first-order form, ${\mathit{q}}^{\prime }=\mathit{u}\left(\mathit{t},\mathit{q}\right)$. For this problem, the system of equations can be rewritten in first-order form using the substitutions $\mathit{w}={\mathit{x}}^{\prime }$ and $\mathit{z}={\mathit{y}}^{\prime }$. With these substitutions, the system contains 28 first-order equations organized into four groups of seven equations:

$$\left[\begin{array}{c}{x}_i^{\prime }\\ y_i^{\prime }\\ w_i^{\prime }\\ z_i^{\prime }\end{array}\right]=\left[\begin{array}{c}{w}_i\\ z_i\\ f_i(x,y)\\ h_i(x,y)\end{array}\right]$$.

The solution vector produced by solving the system has the form

$$\mathit{q}=\left[\begin{array}{c}{\mathit{x}}_{\mathit{i}}\\ {\mathit{y}}_{\mathit{i}}\\ {\mathit{w}}_{\mathit{i}}\\ {\mathit{z}}_{\mathit{i}}\end{array}\right].$$

Therefore, writing the original equations in terms of the solution vector $\mathit{q}$ yields

$$\left[\begin{array}{c}{{\mathit{x}}_{\mathit{i}}}^{\prime }\\ {{\mathit{y}}_{\mathit{i}}}^{\prime }\\ {{\mathit{w}}_{\mathit{i}}}^{\prime }\\ {{\mathit{z}}_{\mathit{i}}}^{\prime }\end{array}\right]=\left[\begin{array}{c}\left({\mathit{q}}_{15},{\mathit{q}}_{16},...,{\mathit{q}}_{21}\right)\\ \left({\mathit{q}}_{22},{\mathit{q}}_{23},...,{\mathit{q}}_{28}\right)\\ {\mathit{f}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)\\ {\mathit{h}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)\end{array}\right],$$

where $\mathit{x}=\left({\mathit{q}}_1,{\mathit{q}}_2,...,{\mathit{q}}_7\right)$ and $\mathit{y}=\left({\mathit{q}}_8,{\mathit{q}}_9,...,{\mathit{q}}_{14}\right)$. With the equations written in the first-order form ${\mathit{q}}^{\prime }=\mathit{u}\left(\mathit{t},\mathit{q}\right)$, you can now write a function that calculates the components during each time step of the solution process. A function that codes this system of equations is:

function dqdt = pleiades(t,q)
  x = q(1:7);
  y = q(8:14);
  xDist = (x - x.');
  yDist = (y - y.');
  r = (xDist.^2+yDist.^2).^(3/2);
  m = (1:7)';
  dqdt = [q(15:28);
          sum(xDist.*m./r,1,'omitnan').';
          sum(yDist.*m./r,1,'omitnan').'];
end

In this function, the $\mathit{x}$ and $\mathit{y}$ values are extracted directly from the solution vector $\mathit{q}$, as are the first 14 components of the output. Then, the function uses the position values to calculate the distances between all seven stars, and these distances are used in the code for ${\mathit{f}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)$ and ${\mathit{h}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)$.

Set Optional Parameters

Use the odeset function to set the value of several optional parameters:

• Specify stringent error tolerances of 1e-13 and 1e-15 for the relative and absolute error tolerances, respectively.

• Turn on the display of solver statistics.

opts = odeset("RelTol",1e-13,"AbsTol",1e-15,"Stats","on");

Solve System of Equations

Create a column vector with the initial conditions and a time vector with regularly spaced points in the range $\left[0,15\right]$. When you specify a time vector with more than two elements, the solver returns solutions at the time points you specify.

init = [3 3 -1 -3 2 -2 2 ...
        3 -3 2 0 0 -4 4 ...
        0 0 0 0 0 1.75 -1.5 ...
        0 0 0 -1.25 1 0 0]';
tspan = linspace(1,15,200);

Now, solve the equations using both ode78 and ode89 by specifying the equations, time span, initial values, and options. Use tic and toc to time each solver for comparison (note that timings can differ depending on your computer).

tic
[t78,q78] = ode78(@pleiades,tspan,init,opts);

14963 successful steps
549 failed attempts
201899 function evaluations

toc

Elapsed time is 2.941335 seconds.

tic
[t89,q89] = ode89(@pleiades,tspan,init,opts);

4181 successful steps
56 failed attempts
68726 function evaluations

toc

Elapsed time is 1.186304 seconds.

The output indicates that ode89 is best suited for solving the problem, because it is faster and has fewer failed steps.

Plot Solution Curves

The first 14 components of q89 contain the x and y positions for each of the seven stars, as obtained by ode89. Plot these solution components to see the trajectories of all the stars over time.

plot(q89(:,1),q89(:,8),'--',...
     q89(:,2),q89(:,9),'--',...
     q89(:,3),q89(:,10),'--',...
     q89(:,4),q89(:,11),'--',...
     q89(:,5),q89(:,12),'--',...
     q89(:,6),q89(:,13),'--',...
     q89(:,7),q89(:,14),'--')
title('Position of Pleiades Stars, Solved by ODE89')
xlabel('X Position')
ylabel('Y Position')

Create Animation of Results

Because the trajectories of the stars overlap considerably, a better way to visualize the results is to create an animation showing the stars move over time. The function AnimateOrbits, included as a local function at the end of this example, accepts the output from a solver for this problem and creates an animated GIF file in the current folder that shows the stars move along their trajectories.

For example, you can generate an animation with the output from the ode89 solver using the command

AnimateOrbits(t89,q89);

The following is a sample output GIF.

References

[1] Hairer, E., et al. Solving Ordinary Differential Equations I: Nonstiff Problems. 2nd rev. ed, Springer, 2009.

[2] “Pleiades.” Wikipedia, 21 June 2021. Wikipedia, https://en.wikipedia.org/wiki/Pleiades.

Local Functions

This section includes the local functions called by the ODE solver to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function dqdt = pleiades(t,q)
  x = q(1:7);
  y = q(8:14);
  xDist = (x - x.');
  yDist = (y - y.');
  r = (xDist.^2+yDist.^2).^(3/2);
  m = (1:7)';
  dqdt = [q(15:28);
          sum(xDist.*m./r,1,'omitnan').';
          sum(yDist.*m./r,1,'omitnan').'];
end
%-----------------------------------------------------------------
function AnimateOrbits(t,q)
for k = 1:length(t)
    plot(q(:,1),q(:,8),'--',q(:,2),q(:,9),'--',...
         q(:,3),q(:,10),'--',q(:,4),q(:,11),'--',...
         q(:,5),q(:,12),'--',q(:,6),q(:,13),'--',...
         q(:,7),q(:,14),'--')
    hold on
    xlim([-20 20])
    ylim([-10 10])
    sz = 15;
    plot(q(k,1),q(k,8),'o','MarkerSize',sz,'MarkerFaceColor','r')
    plot(q(k,2),q(k,9),'o','MarkerSize',sz,'MarkerFaceColor','k')
    plot(q(k,3),q(k,10),'o','MarkerSize',sz,'MarkerFaceColor','b')
    plot(q(k,4),q(k,11),'o','MarkerSize',sz,'MarkerFaceColor','m')
    plot(q(k,5),q(k,12),'o','MarkerSize',sz,'MarkerFaceColor','c')
    plot(q(k,6),q(k,13),'o','MarkerSize',sz,'MarkerFaceColor','y')
    plot(q(k,7),q(k,14),'o','MarkerSize',sz,'MarkerFaceColor',[210 120 0]./255)
    hold off
    drawnow
    M(k) = getframe(gca);
    im{k} = frame2im(M(k));
end

filename = "orbits.gif";
for idx = 1:length(im)
    [A,map] = rgb2ind(im{idx},256);
    if idx == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
    end
end
close all
end