Trouble with ODE45, (not enough data of a planet's orbit)
2 views (last 30 days)
Show older comments
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
Accepted Answer
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
Jan
on 18 Dec 2016
Edited: 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. ;-)
More Answers (2)
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
Fabio Ricci
on 19 Dec 2016
1 Comment
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.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!