MATLAB Answers

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

12 views (last 30 days)
Fabio Ricci
Fabio Ricci on 18 Dec 2016
Commented: Jan on 19 Dec 2016
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

  0 Comments

Sign in to comment.

Accepted Answer

Mischa Kim
Mischa Kim on 18 Dec 2016
Edited: Mischa Kim on 18 Dec 2016
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.

  7 Comments

Show 4 older comments
Image Analyst
Image Analyst on 18 Dec 2016
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
Fabio Ricci on 18 Dec 2016
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
Jan on 18 Dec 2016
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. ;-)

Sign in to comment.

More Answers (2)

Jan
Jan on 18 Dec 2016
Edited: Jan on 18 Dec 2016
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]

  0 Comments

Sign in to comment.


Fabio Ricci
Fabio Ricci on 19 Dec 2016
OK you are right, it's my fault thank you for the help.

  1 Comment

Jan
Jan on 19 Dec 2016
@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.

Sign in to comment.

Sign in to answer this question.


Translated by