TwoBodyDemo.m

Contents

Overview

This script illustrates the use of ode23tx for the solution of the two-body problem

$$ x''(t) = -\frac{x(t)}{r(t)^3}, $$

$$ y''(t) = -\frac{y(t)}{r(t)^3}, $$

where

$$ r(t) = \sqrt{x(t)^2 + y(t)^2}. $$

Here $x$ and $y$ are the coordinates of a small body (such as a space ship) relative to a (fixed) large body (such as a planet), and $t$ is time. In other words, the parametric curve $C: t \mapsto (x(t),y(t))$ describes the orbit of the small body around the large body.

In order to apply an IVP solver we need to convert the system of two second-order ODEs to a system of four first-order ODEs, i.e.,

$$ \mathbf{y} = \left[y_1, y_2, y_3, y_4 \right]^T = \left[x, x', y, y' \right]^T $$

so that

$$ \mathbf{y}' = \left[y'_1, y'_2, y'_3, y'_4 \right]^T = \left[y_2, -\frac{y_1}{r^3}, y_4, -\frac{y_3}{r^3} \right]^T $$

This first-order system is coded in the function twobody.m.

Beginning of code

Initialize

clear all
close all

Example 1

Use ode23tx to solve the two-body problem with initial condition

$$ x(0) = 1,\ x'(0) = 0,\ y(0) = 0,\ y'(0) = 1, $$

i.e., with the small body starting at $(1,0)$ with unit initial velocity in the "upward" direction. The large body is to be interpreted as a point mass located at the origin. We use a time interval of $[0,2\pi]$. The horizontal axis in the plot is the $t$-axis. The different colors correspond to

ode23tx(@twobody,[0 2*pi],[1; 0; 0; 1]);
pause

To get a phase plane plot, i.e., a plot of the $t$-parametrized actual orbit determined by the first and third components of the system of four first-order ODEs we store the output of ode23tx in t and y and then plot those first and third components of y together.

y0 = [1; 0; 0; 1];
[t,y] = ode23tx(@twobody,[0 2*pi],y0);
plot(y(:,1),y(:,3),'-',0,0,'ro')
xlabel('x'), ylabel('y')
axis equal
pause

Example 2

Now we use ode23tx to solve the two-body problem with slightly modified initial conditions

$$ x(0) = 1,\ x'(0) = 0,\ y(0) = 0,\ y'(0) = 1.3, $$

i.e., we increase the initial velocity in the "upward" direction. We also change the time interval to $[0,12\pi]$. Again, the horizontal axis in the plot is the $t$-axis, and the different colors correspond to

y0 = [1; 0; 0; 1.3];
ode23tx(@twobody,[0 12*pi],y0);
pause

The corresponding phase plane plot for this modified $y$-velocity is given by

y0 = [1; 0; 0; 1.3];
[t,y] = ode23tx(@twobody,[0 12*pi],y0);
plot(y(:,1),y(:,3),'-',0,0,'ro')
xlabel('x'), ylabel('y')
axis equal


© Greg Fasshauer, IIT. Published with MATLAB® 7.12