Trouble with ODE45, (not enough data of a planet's orbit)
2 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
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 件のコメント
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
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]
0 件のコメント
Fabio Ricci
2016 年 12 月 19 日
1 件のコメント
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!