Orbit by euler's method

Hello, I need to create a script that uses these iteration functions to create an orbit chart, but all the way I tried the most I could get was straight, thanks for the help.

Here is one of the codes I tried to do.

clc
clear
v(1,1)=0;
v(1,2)=1;
r(1,1)=20;
r(1,2)=0;

i=2;
delt=0.01;
t=0
while t<=200
    mr=sqrt(((r((i-1),1))^2)+((r((i-1),2))^2));
    v(i,:)=v((i-1),:)+((-r((i-1),:))/(mr^3))*delt;
    r(i,:)=r((i-1),:)+(((v((i-1),:)+v(i,:))/2)*delt);
    t=t+delt
    i=i+1;
end
plot(r(:,1),r(:,2));

Answer 1

[James Tursa]

I don't see anything obviously wrong with your code. So other possibilities to consider are:

- The step size is too large for the dynamics involved

- The initial velocity is too large and you are in a hyperbolic orbit (looks almost linear)

- You are missing a constant factor on the acceleration term (the "a" part) of the velocity equation

For instance, run this version of your code (using the same equations) with a smaller initial velocity:

clc
clear
delt=0.01;
t=0;
tend = 4000;
N = floor(tend/delt) + 1;
r = zeros(N,2);
v = zeros(N,2);

v(1,1)=0;
v(1,2)=1;
r(1,1)=20;
r(1,2)=0;

v(1,2)=0.1; % Force a smaller initial velocity

i=2;

while t<=tend
    mr=sqrt(((r((i-1),1))^2)+((r((i-1),2))^2));
    v(i,:)=v((i-1),:)+((-r((i-1),:))/(mr^3))*delt;
    r(i,:)=r((i-1),:)+(((v((i-1),:)+v(i,:))/2)*delt);
    t=t+delt;
    i=i+1;
end
plot(r(:,1),r(:,2));
maxr = max(r);
minr = min(r);
aver = (minr + maxr)/2;
d = max(maxr - minr)/2;
xlim([aver(1)-d aver(1)+d]);
ylim([aver(2)-d aver(2)+d]);
axis square
grid on

And you get this, which looks more like what you were probably expecting:

In this case you can see that the Euler scheme has a systematic integration error that puts the trajectory too high at each step, resulting in an outward spiral instead of the correct closed loop ellipse.

[Lucas Diniz]

Thank you very much. The teacher must have been wrong about the initial values. :)

[James Tursa]

I would definitely double check that acceleration term. Currently you just have a 1 in the numerator, but normally I would have expected to see a constant factor there representing the attracting force constant (unless maybe you are using some type of normalized units).

[Lucas Diniz]

I have a PDF with the deduction of these functions, but it is in portuguese. If you are interested, start on page 13.

[Lucas Diniz]

Again, thanks for the help.

[James Tursa]

Assuming some type of normalized units resulting in that 1 in the numerator of acceleration, the escape velocity of your system at radius r is sqrt(2/r). So for your first case where you start at r=20, the escape velocity is sqrt(2/20) = 0.316. Your initial velocity of 1 for that case is well above escape velocity, so it is in a hyperbolic orbit, hence the somewhat linear looking graph when it is plotted instead of an ellipse. Looks like you're going to have problems with those other cases as well.

[Efan Nonti]

James Tursa can i asking? What is the v(i,:) & r(i,:) mean?

[James Tursa]

r(i,:) is the 2D position vector at the i'th time, and v(i,:) is the 2D velocity vector at the i'th time.