Problem with a little solar system simulation

5 views (last 30 days)
Jorge
Jorge on 30 Jun 2011
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. :)

Answers (0)

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!