Issues solving an ODE system for the n-body problem (Edited)
이전 댓글 표시
I need to solve the ODE system of the n-body problem (in 2D) in order to get future positions and velocities from given initial values of the system and afterwards make a simulation of their behavior.
The ODE system is

where r are the positions, v the velocities, m the masses and Fg(r) the gravitational force.
I use ode45(odefun,tspan,y0). I need to define a function, I call it ode_nbody, which holds the mentioned system and is used in this fashion ode45(ode_nbody,tspan,y0). ode_nbody receive as a parameter the masses. The function I've written is shown below
function [ddt] = ode_nbody(t, y, m)
% y = [r1x, r2x, ..., rnx, r1y, r2y, ..., rny, v1x, v2x, ..., vnx, v1y, v2y, ..., vny]
G = 6.67408E-11;
n = length(m);
r = [y(1:n), y(n+1:2*n)]'; % Vector of positions that makes easier the acceleration computation
for i = 1:n
acc(:,i) = zeros(2,1);
for j = 1:n
if i ~= j
acc(:,i) = acc(:,i) + ...
G * (m(j) .* (r(:,j) - r(:,i)))./(norm(r(:,j) - r(:,i))^3);
end
end
end
% ddt = [v1x, v2x, ..., vnx, v1y, v2y, ..., vny, a1x, a2x, ..., anx]
ddt(1:2*n) = y(2*n+1:4*n);
ddt(2*n+1:4*n) = [acc(1,1:n) acc(2,1:n)];
ddt = ddt';
end
The results is

As you may have noticed, all the planets converge to the sun (that yellow "o" in the center) which, obviusly was not expected.
댓글 수: 4
I recovered the removed content from the Google cache (something which anyone can do). Editing away your question is very rude. Someone spent time reading your question, understanding your issue, figuring out the solution, and writing an answer. Now you repay that kindness by ensuring that the next person with a similar question can't benefit from this answer.
M. Ezequiel
2021년 7월 20일
Rik
2021년 7월 20일
Once an answer has been posted, you are no longer able to delete a question. I have understood Mathworks has (or will soon) added a remark explaining this when you post a question. If you did not intend to be rude, then there is no longer an issue, simply leave the question as-is.
If you have already solved it, feel free to post the solution as an answer and accept it. That way you might help someone in the future with a similar issue.
M. Ezequiel
2021년 7월 20일
답변 (1개)
Walter Roberson
2021년 7월 18일
0 개 추천
You did not say what the difficulty is ?
Note that you are only tracking the orbits for 10 seconds
댓글 수: 1
Rik
2021년 7월 20일
Oh yep, sorry! The difficulty is that if you put a bigger interval, for example [0 1000000], the trajectories converge to the sun (the yellow "o" in the image) but they're suppose to depict a circular orbit.

카테고리
도움말 센터 및 File Exchange에서 Earth and Planetary Science에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!