Matlab simulation for planet motion

10 次查看（过去 30 天）
Niklas Kurz
Niklas Kurz 2022-3-16
评论： Niklas Kurz 2022-3-17
There were some attemps simulating planetary motion already, but I think mine is straightforward by solving and updating position via with Euler Cromers method:
t = 0;
while t < 10
pos1 = [1 2 3];
pos2 = [4 5 6];
m1 = 1;
m2 = 2;
G = 1;
r1 = pos1-pos2;
r2 = pos2-pos1;
F1 = G*m1*m2/norm(r1).^2.*r1/norm(r1);
F2 = G*m1*m2/norm(r2).^2.*r2/norm(r2);
dt = 0.1;
p1 = [0 100 0];
p2 = [0 100 0];
p1 = p1+F1.*dt;
p2 = p2+F2.*dt;
pos1 = pos1+p1/m1;
pos2 = pos2+p2/m2;
t = t+dt;
hold all;
plot3(pos1(1),pos1(2),pos1(3),'rx')
plot3(pos2(1),pos2(2),pos2(3),'bx')
end
However I don't really receive a plot of multiple data points, just 2 crosses remaining stationary. Also I get a 2-D plot even though I reverted to plot3
  1 个评论
KSSV
KSSV 2022-3-16
You can change it to 3D using view.
plot3(pos1(1),pos1(2),pos1(3),'rx')
plot3(pos2(1),pos2(2),pos2(3),'bx')
view(3)

请先登录，再进行评论。

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

采纳的回答

James Tursa
James Tursa 2022-3-16
The initial condition for position and velocity need to be outside the loop, prior to loop entry.

更多回答（1 个）

KSSV
KSSV 2022-3-16
Ran in:
t = 0;
m1 = 1;
m2 = 2;
G = 1;
pos01 = [1 2 3];
pos02 = [4 5 6];
pos1 = zeros([],3) ;
pos2 = zeros([],3) ;
iter = 0 ;
while t < 10
iter = iter+1 ;
r1 = pos01-pos02;
r2 = pos02-pos01;
F1 = G*m1*m2/norm(r1).^2.*r1/norm(r1);
F2 = G*m1*m2/norm(r2).^2.*r2/norm(r2);
dt = 0.1;
p1 = [0 100 0];
p2 = [0 100 0];
p1 = p1+F1.*dt;
p2 = p2+F2.*dt;
pos1(iter,:) = pos01+p1/m1;
pos2(iter,:) = pos02+p2/m2;
pos01 = pos1(iter,:) ;
pos02 = pos2(iter,:) ;
t = t+dt;
end
figure
hold on
plot3(pos1(:,1),pos1(:,2),pos1(:,3),'rx')
plot3(pos2(:,1),pos2(:,2),pos2(:,3),'bx')
view(3)
  1 个评论
Niklas Kurz
Niklas Kurz 2022-3-17
This looks much better to me regarding number of points. Nevertheless there is still something weird with the coding going around. Even if it should equal the code that I assimilated from Glowscript:
G = 1
star = sphere(pos = vector(0,0,0), radius = 0.2, color = color.yellow, mass = 1000, momentum = vector(0,0 ,0), make_trail=True)
plan = sphere(pos = vector(1,0,0), radius = 0.5, color = color.blue , mass = 1 , momentum = vector(0,30,0), make_trail=True)
while (True):
rate(500)
r_star = star.pos - plan.pos
r_plan = plan.pos - star.pos
star.force = -G*star.mass*plan.mass/(mag(r_star)**2)*(r_star)/mag(r_star)
plan.force = -G*star.mass*plan.mass/(mag(r_plan)**2)*(r_plan)/mag(r_plan)
star.momentum = star.momentum + star.force * dt
plan.momentum = plan.momentum + plan.force * dt
star.pos = star.pos + star.momentum/star.mass * dt
plan.pos = plan.pos + plan.momentum/plan.mass * dt
t = t+dt

请先登录，再进行评论。

类别

在 Help Center 和 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!

Translated by