Problem with a little solar system simulation
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Earth and Planetary Science en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!