How can i make my data more efficient in time when running script?

I'm calculating the orbit of Earth and other planets around the sun using force of gravity. The script i have runs quick enough using a single planet, but it then begins to take longer as i add more planets (obviously).

I'm hoping if someone can give me some help on how to use the zeros function or another way to make the time taken to process the calculations more efficiently. I'm using a for loop to calculate and update each position through each iteration.

here is an example of the data:

days=365*48; % number of steps to make  
stepsize=1/48; %(dt) or in Days  
sunm =1.9891e30; % Sun Mass (kg)  
G=1.4879e-34; % AU^3/(kg*Day^2)  
%%Earth Initial Conditions  
Ex=-9.320752400360321E-01; % Earth pos at x AU  
Ey=-3.683263876537671E-01;  % Earth pos at y AU  
Ez=1.235087446387007E-05; % Earth pos at z AU  
E_Vx=6.036562219130709E-03;% Earth x velocity AU/day  
E_Vy=-1.606778488778610E-02;% Earth y velocity AU/day  
E_Vz=6.657490105653882E-07; % Earth z velocity AU/day

%%Earth orbit  
for count=1:days      
    % Save the current position of earth into a matrix
    finalposition(count,1)=Ex;
    finalposition(count,2)=Ey;
    finalposition(count,3)=Ez;

      %Calculate current Acceleration
      R=norm([Ex,Ey,Ez]); 

      % assuming the sun is always at [0,0,0]
      Atot=G*sunm/R^2;
      unitvector=[Ex,Ey,Ez]/R;

      Ax=-Atot*unitvector(1); %acceleration in the x
      Ay=-Atot*unitvector(2); %acceleration in the y
      Az=-Atot*unitvector(3); %acceleration in the z

      %Update the position of Earth
      Ex=Ex+E_Vx*stepsize+0.5*Ax*stepsize^2;
      Ey=Ey+E_Vy*stepsize+0.5*Ax*stepsize^2;
      Ez=Ez+E_Vz*stepsize+0.5*Ax*stepsize^2;

      %Update the velocity of Earth
      E_Vx=E_Vx+Ax*stepsize;
      E_Vy=E_Vy+Ay*stepsize;
      E_Vz=E_Vz+Az*stepsize;
  end

hold on
plot3(finalposition(:,1),finalposition(:,2),finalposition(:,3),'b')

Thank you very much for your time, Alex :)

[Thorsten]

How are R and sunm defined?

[Alex]

Sorry, edited. I'm assuming sun is a single point at [0, 0, 0]. sunm (sun mass) is 1.9891e30 (kg). G is the gravitational constant.

[Titus Edelhofer]

Hi,

I would strongly suggest to use a "real" ODE solver instead of the update you are using. Computing the orbit requires quite some good precision to give meaningful results. Take a look at e.g. ode45 ...

Titus

Answer 1

[the cyclist]

Most importantly, not for speed, but for correctness, is that it looks like your formulas for "Update the position of the Earth" incorrectly use Ax for the y and z components.

Probably the biggest thing to help the speed will be to preallocate finalposition. Put this line ahead of your for loop:

finalposition = zeros(days,3);

This assigns memory for that variable ahead, rather than growing that array.

For speed, you could also just define one vector A, and do your updates like

A=-Atot*unitvector

Ditto for the other components.

[Alex]

Thank you for picking up that error, and thank you for your suggestions.