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

1 次查看（过去 30 天）
Fabio Ricci
Fabio Ricci 2016-12-18
评论： Jan 2016-12-19
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

请先登录，再回答此问题。

采纳的回答

Mischa Kim
Mischa Kim 2016-12-18
编辑：Mischa Kim 2016-12-18
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 个评论
显示 5更早的评论
Fabio Ricci
Fabio Ricci 2016-12-18
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 2016-12-18
编辑：Jan 2016-12-18
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. ;-)

请先登录，再进行评论。

更多回答（2 个）

Jan
Jan 2016-12-18
编辑：Jan 2016-12-18
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]
Fabio Ricci
Fabio Ricci 2016-12-19
OK you are right, it's my fault thank you for the help.
  1 个评论
Jan
Jan 2016-12-19
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by