Problem with a little solar system simulation
5 views (last 30 days)
Show older comments
HI!
i'm starting using Matlab, and, as a little exercise im doing a simulation that only includes the Sun and the Earth, so my problem is the next: im using the Newton 2º Law to solve it and get:
r'' = G*M/r^2 - r'^2/r
where r is the distance Sun-Earth, M is the mass of the sun and G the gravitational constant.
In order to solve that, i have written two scripts:
clear all
G=6.67428*10^(-11); %m³/(kg·s²)
Msun=1.9891*10^30; %kg
t_0=0; %s initial time
t_f=63936000; %s final time
N=15000; %steps
M_t=5.98*10^24; %kg
v0_t=29440; %m/s
r0_t=149600000000; %m
omega0_t=v0_t/r0_t; %initial angular velocity
theta0_t=0; %inicial angle
g=@(r_t,v_t,theta_t,omega_t,t_t)(v_t);
h=@(r_t,v_t,theta_t,omega_t,t_t)(G*Msun./(r_t.^2) - r_t.*omega_t.^2);
j=@(r_t,v_t,theta_t,omega_t,t_t)(omega_t);
k=@(r_t,v_t,theta_t,omega_t,t_t)(2./r_t.*v_t.*omega_t);
[r_t,v_t,theta_t,omega_t,t_t]=int_rk2b(g,h,j,k,r0_t,v0_t,theta0_t,omega0_t,t_0,t_f,N);
plot(t_t,r_t)
the second script:
function[r,v,theta,omega,t]=int_rk2b(g,h,j,k,r0,v0,theta0,omega0,t0,tf,N)
dt=(tf-t0)/(N-1);
t=linspace(t0,tf,N);
r=zeros(1,N);
v=zeros(1,N);
theta=zeros(1,N);
omega=zeros(1,N);
t(1)=t0;
r(1)=r0;
v(1)=v0;
theta(1)=theta0;
omega(1)=omega0;
for i=1:N-1
tmid=t(i)+dt/2;
rmid=r(i)+g(r(i),v(i),theta(i),omega(i),t(i))*dt/2;
vmid=v(i)+h(r(i),v(i),theta(i),omega(i),t(i))*dt/2;
thetamid=theta(i)+j(r(i),v(i),theta(i),omega(i),t(i))*dt/2;
omegamid=omega(i)+k(r(i),v(i),theta(i),omega(i),t(i))*dt/2;
r(i+1)=r(i)+g(rmid,vmid,thetamid,omegamid,tmid)*dt;
v(i+1)=v(i)+h(rmid,vmid,thetamid,omegamid,tmid)*dt;
theta(i+1)=theta(i)+j(rmid,vmid,thetamid,omegamid,tmid)*dt;
omega(i+1)=omega(i)+k(rmid,vmid,thetamid,omegamid,tmid)*dt;
end
And that's all, with the first script i call the second one to solve these second order differential equations but if you run the program, results are completly wrong i dont know why, i hope anyone can find out my mistake. :)
0 Comments
Answers (0)
See Also
Categories
Find more on Earth and Planetary Science 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!