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.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!