Trouble with ODE45, (not enough data of a planet's orbit)

Hi everybody, i was simulating the orbit of Uranus with Matlab with the initial condition taken from there: http://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html

With this following code: Function:

function [ Ydot ] = f(t,Y)
global m1 G %Masse dei tre corpi
x0=Y(1:3);
d0=(x0)/norm(x0)^3;
Ydot(1:3)=Y(4:6);
Ydot(4:6)=-G*m1*d0;
Ydot=Ydot(:);
end

Body:

global m1 G %Masse corpi 

m1=2e30;
G=6.67e-11;
Perkm=2.741e12; %m
Velperi=7.11e3;%m/s
x00=[Perkm;0;0]; xp0=[0;Velperi;0];

options=odeset('RelTol',1e-10,'AbsTol',1e-10);
[T1,Y1]=ode45('fUranus',[0 3e+12],[x00;xp0],options);

The problem is the output of Y1, which is in cartesian coordinates (X,Y,Z) where Z is alwais 0. I expect that there is a quite gradual decrease of the X, because of the elliptical orbit but i obtain something like a Jump

but when i try to plot the orbit i obtain a complite ellipse

Thank you in advance. Fabio Ricci

Answer 1

[Mischa Kim]

Fabio, you control the number of data points that are returned in

 [T1,Y1] = ode45('fUranus',tspan,[x00;xp0],options);

by specifying tspan as

 tspan = linspace(0, 3e+12, N);

where N is the number of data points.

Note that this does not impact the accuracy of our results, which is set in odeset.

[Fabio Ricci]

Thank you a lot, I was missing this information.

[Fabio Ricci]

Sorry I belived that the problem was solved, now I have more data points but the jump persists.

[Fabio Ricci]

I've noticed that when i use:

tspan=linspace(0,3e+12,3000);

I obtain a vector that is not equally spaced

[Jan]

Do not post tiny screenshots. Posting data as text is easier to read. Why do you assume, that linspace does not reply equally spaced values? Of course, there is a little rounding, but due to the limited precision of floating point values, this cannot be avoided.

tspan = linspace(0,3e+12,3000)
format long g
unique(diff(tspan))

Looks fine for me.

I do not see, what the problem in the original question is. I do not see, what you call "jump".

[Image Analyst]

When my mouse hovers over the screenshot, Firefox blows it up to full screen. Not sure why it doesn't on yours. Maybe I've got some plugin you don't.

[Fabio Ricci]

Thank you a lot for the help; With Chrome evrything looks perfect, atleast for me, by the way i am going to upload both the data file both the image to avoid problems.

Now, the problem is the non equaly spaced vector created with linspace, sounds like a contraddiction.

i made a try with

tspan=linspace(0,10,11);

and everything is fine.

but when i tried with

tspan=linspace(0,3e+12,13);

i obtained

As you can see from the data or the image the data points are not equally spaced, in fact they start from e+11.

[Jan]

The posted values are equally distributes with a step size or 2.5e11. Perfect. You divide the interval from 0 to 3e12 in 13 points, which means 12 intervals, then the step size is 3e12/12, which gives 2.5e11. The first point is 0, the second is 2.5e11, the third is 2*2.5e11=5e11 and so on. Everything is fine.

What do you expect instead???

Again: The data in the original question do not only look nice, even the created diagram shows, that there are no jumps. Perhaps you are confused by the "eXYZ" notation?

I have to click on the screenshot in my Firefox/Windows7, but then the current page disappears. I do not have an "open the figure in a new tab" field in the context menu. Therefore inspecting the screenshot and reading the post is tedious. Showing the OP, that the data are smooth would be easy, if I could copy&paster the data.

[EDITED] I've found a solution: The context menu offers to set the image as desktop background. Now I can see it on my second monitor while the thread is still open. ;-)

Answer 2

[Jan]

To post is as answer also: The shown data do not contain any jumps. Everything is fine. (It would be easier to show this, if you had posted the data as text and not as screenshot - copy&paste is the simplest way to work with data from the forum.)

x = [2.5867e10, 1.4022e10, 2.177e9, -9.6684e9, -1.9635e10, -2.9602e10]
format shortG
diff(x)
>> -1.1845e+10  -1.1845e+10  -1.1845e+10  -9.9666e+09   -9.967e+09

Do you see? The steps between the data are in the same magnitude near to -1e10. No jumps. There is no problem and therefore no solution.

If you still do not see it, divide the data by 1e10:

x = [2.5867e10, 1.4022e10, 2.177e9, -9.6684e9, -1.9635e10, -2.9602e10] / 1e10
>> [2.5867   1.4022    0.2177   -0.96684   -1.9635   -2.9602]

Answer 3

[Fabio Ricci]

OK you are right, it's my fault thank you for the help.

[Jan]

@Fabio: You are welcome.

You decided that Mischa's answer has solved your problem. Perhaps I am the one who did not understand, what the problem is. Sigh.