Get x and y coordinates of motion from ode45

Hi,

I would like to plot the motion of a satellite in the x-y-plane around the Earth. For this motion I have the differential equation

. As initial conditions, I have and with two numerical values for r and v and I want to pass them to ode45 via

[t1, y1]=ode45('unperturbed', [0 100], [r0 v0]); 

How does unperturbed have to look like so it can handle the vectors and returns a position vector for every time t?

So far I have it like this which seems to return valid y coordinates but I don't know how to get the x part.

function dy = unperturbed(t,y)
M=5.594e24;
G=6.674e-11;
rE=6378.137;
r0=rE+800;

dy=zeros(2,1);
dy(1) = y(2);
dy(2) = -G*M*(dy(1)/r0^3);
end

I hope you can help me with this

[Torsten]

How does your equation from above look if you write it out ?

Is it

d^2(r_x)/dt^2 = 0
d^2(r_y)/dt^2 = -g*M*r_y/r0^3

with initial conditions

r_x(0) = r0, d(r_x)(0)/dt = 0, r_y(0) = 0, d(r_y)(0)/dt = v0

?

Answer 1

[James Tursa]

You've got a 4th order system (2 equations for x and y, each 2nd order --> 2 x 2 = 4), so your state vector will need to be 4 elements, not 2 elements like you currently have. Define the state vector as follows and then work from there:

y(1) = the x position
y(2) = the y position
y(3) = the x velocity
y(4) = the y velocity

Then in your derivative function you will have something like this:

dy=zeros(4,1); % 4th order system so 4 elements
dy(1) = y(3);
dy(2) = y(4);
% the -GM *rvec/r^3 calculation goes here (this will be a vector result)
dy(3) = _____; % you fill this in, based on the vector result above
dy(4) = _____; % you fill this in, based on the vector result above

Side Note: The DEQ in your Question is missing a minus sign (even though you do have it in your code)

[Backtobasics]

Thank you very much, I think I understand now. I'll try to get it running later and I'll let you know. And there was a minus sign missing indeed, I corrected it.

[Backtobasics]

Okay, maybe I was too fast to proclaim "I understand" :D

This is my code in the main file now

% initial conditions
r=rE+800;
v=sqrt(2*mu/r);

y(1)=r;
y(2)=0;
y(3)=0;
y(4)=v;

% part A

[t1, y1]=ode45('unperturbed', [0 100], y);
plot(y1(1,:), y1(2,:));

and in the fuction file it's

function dy = unperturbed(t,y)
dy=zeros(4,1);
dy(1) = y(3);
dy(2) = y(4);

% the -GM *rvec/r^3 calculation goes here (this will be a vector result)
M=5.594e24;
G=6.674e-11;
rE=6378.137;
r0=rE+800;
vec=-G*M*[r0; 0]/(norm([r0; 0])^3); %dummy, I didn't know how to call it instead

dy(3) = vec(1);
dy(4) = vec(2);
end

This leads to results that don't really work. What am I doing wrong? I've got to say that this my first time working with ode's so I'm still trying to really understand how they work.

[James Tursa]

The rvec and r in the gravity calculation should be based on the incoming state vector, not on fixed values as you have done. E.g.,

vec = -G*M*y(1:2)/(norm(y(1:2))^3); % <-- gravity force is based on current state vector position

And actually, you don't really even need the vec temporary variable. E.g.,

dy(1:2) = y(3:4);
dy(3:4) = -G*M*y(1:2)/(norm(y(1:2))^3);

[Backtobasics]

Okay, thank you very much! So my idea wasn't completely wrong.

Plotting plot(y1(1,:), y1(2,:)) should plot the x and y coordinates and should result in a circular shape, right? Because what I get is the following:

[James Tursa]

The expected motion is dependent on initial velocity. You have not started with a circular velocity. E.g.,

v=sqrt(2*mu/r); % <-- Escape velocity
v=sqrt(mu/r); % <-- Circular velocity

To see a circle on your plot you will probably want to set the x and y axis limits the same, and use axis square.

[Backtobasics]

True, I used the wrong formula. And I also realized that sometimes I used km and sometimes m so I thought this might have been the problem - turns out it is not. Using your axis properties, it looks like this now.

So how I read this is that the objects seems to start at r=[7178; 7178] (while it should be at r=[7178; 0]) and then falls straight into the "center of gravity". But why?

Sorry to bother you with this, probably it's a tiny fault somewhere but I went through the code a couple of times now and I cannot find it. For completeness I'll attach the complete files to this, maybe this helps.

[James Tursa]

This is what you have in your current code (1st & 4th):

plot(y1(1,:), y1(4,:));

But this is what you had in your earlier code (1st & 2nd):

plot(y1(1,:), y1(2,:));

Why the change? The x and y position values are in the first two elements of the state vector.

Also, for consistency, note these lines:

M=5.594e24;
G=6.674e-11;
mu=3.986e14;

I would either use M and G throughout your code, or mu throughout your code, but I would not use both since one is dependent on the other. E.g., if you want to use M and G, then do this:

mu = G*M;

In fact, in all problems with units, I would strongly advise listing the units and definition on the very line where the constant is defined so that you can easily spot unit mismatches and know what you were doing months later when you look at this code. E.g.,

km2m = 1000; % kilometers to meters
M = 5.594e24; % (kg), mass of earth
G = 6.674e-11; % (m^3/kg*s^2), gravitational constant
mu = G*M; % (m^3/s^2), mu of Earth
rE = 6378.137 * km2m; % (m), equatorial radius of Earth
% initial conditions
altitude = 800 * km2m; % (m), altitude of orbit
r = rE + altitude; % (m), radius of orbit

Another change I would make is to pass in any parameters into your derivative function instead of having multiple definitions of them scattered throughout your code. E.g., pass M and G into your derivative function like this:

function dy = unperturbed(t,y,M,G)
dy = zeros(4,1);
dy(1:2) = y(3:4);
dy(3:4) = -G*M*y(1:2)/(norm(y(1:2))^3);
end

And then call it this way:

[t1, y1] = ode45(@(t,y)unperturbed(t,y,M,G), [0 10000], y);

[Backtobasics]

Oh yes, sorry, I uploaded the wrong file. In this, I tried to combine different y1( )'s to see if maybe y1(1,:) and y1(2,:) are not the ones I have to use. And I definitely see the point in the units problems, thank you.

I found the mistake now and it's silly: Instead of plot(y1(1,:), y1(2,:)), it's plot(y1(:,1), y1(:,2)).

Now it looks a lot better, the only thing left is to find out why the satellite is spiralling down instead of a circular orbit.

[James Tursa]

It's likely spiraling down because of ode45 integration errors. I.e., the ode45 method isn't particularly well suited for this problem. Try tightening the tolerances (see the options structure argument).

[Backtobasics]

Amazing, it worked perfectly! Thank you so, so much!