CODEWITHSUNDEEP
matlabodedifferential-equationsnumerical-integration
Ridgy
Ridgy

Reputation: 13

Basic 2-body interaction using Matlab's ode45

I'm trying to model basic gravitation for an object of negligible mass around a massive body. I've followed the examples provided in the ODE suite documentation, but the results I'm getting are plainly ridiculous.

Here's the function I'm calling with ode45:

function dy = rigid(t,y)
dy = zeros(4,1);

%Position
xx=y(1);
yy=y(2);

%Radius
r=(xx.^2+yy.^2).^0.5;

%Constants
M=10^30;
G=6.67*10^-11;

%dX/dt
dy(1)=y(3); %vx
dy(3)=-M.*G.*xx.*r.^-3; %ax

%dY/dt
dy(2)=y(4); %vy
dy(4)=-M.*G.*yy.*r.^-3; %ay

And here are the solver lines: 

%Init
x_0=1;
y_0=1;
vx_0=0;
vy_0=5;

%ODEs
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);

%Vectors
x=Y(:,1);
y=Y(:,2);

%Plot
plot(x,y)
xlabel('x');
ylabel('y');
title('y=f(x)');

I get a linear plot at the end. Even with initial speed, the position doesn't change over a period of several steps. The only thing I can think of is that I've misunderstood the way to set up my system of ODEs.

I've been pondering this for a while now, and I'm really short on ideas having had done a few searches on the web.

Upvotes: 1

Views: 6317

Answers (1)

A. Donda
A. Donda

Reputation: 8476

Summary: There are problems in integrating Hamiltonian systems with normal numerical integrators, and your special initial conditions aggravate this to the point where the numerical solution has no resemblance with the correct one.

There's nothing wrong with your implementation per se. However, the initial conditions you use are not the best. The constants G and M you use are in SI units, which means the coordinates are in m, the speeds are in m/s, and time is in s. These lines

x_0=1;
y_0=1;
vx_0=0;
vy_0=5;
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);

therefore mean that you are asking for an orbit with a radius of about 1.4 meters and an orbital speed of 5 m/s, and you want this orbit over a period of 17 minutes. Imagine there actually was an object just meters away from a mass of 10^30 kilograms!

So let's try to ask for something more realistic, similar to Earths' orbit, and look at it over 1 year:

x_0=149.513e9;
y_0=0;
vx_0=0;
vy_0=29.78e3;
[T,Y] = ode45(@rigid,[0 31.536e6],[x_0 y_0 vx_0 vy_0]);

And the result is

which looks as expected.

But there is a second problem here. Let's look at this orbit over a period of 10 years ([0 315.36e6]):

Now we no longer get a closed orbit, but a spiral! That's because the numerical integration proceeds with limited precision, and for this set of equations this leads (physically speaking) to a loss of energy. The precision can be increased using parameters to ode45, but ultimately the problem will persist.

Now let's go back to your original parameters and have a look at the result:

This "orbit" is a straight line towards the origin (the sun). Which could be ok, since a straight oscillation is a possible special case of an elliptic orbit. But looking at the coordinates over time:

We see that there is no oscillation. The "planet" falls into the sun and stays there. What's happening here is the same effect as with the larger orbit: Imprecise integration leads to a loss of energy. Moreover, the numerical integration gets stuck: We asked for a period of 1000 s, but the integration does not proceed beyond 1.6e-10 seconds.

As far as I can tell, neither ode45 nor any of the other standard Matlab integrators are adequate to deal with this problem. There are special numerical integration algorithms designed to do so, specifically for Hamiltonian systems, called symplectic integrators. There is a file exchange entry that provides different implementations. Also see this question and answers for more pointers.

Upvotes: 2

Related Questions