problem ode45 system of differential equation

hi, i have to solve this problem about the orbit of the halley comet, but when i run the code i get unphysical values.

up on this text , i uploaded the two equation to solve, and C= -M*G.

i splitted the two equation of the acceleration in four eqs.

anyone can tell me what's the problem with the code?

CODE:

tSpan=[0 2.33e9]

y0=[8.6e7 0 0 55.1]

[tSol,ySol]=ode45(@cometa,tSpan,y0)

plot(tSol,ySol(:,1));

hold on

plot(tSol,ySol(:,2));

plot(tSol,ySol(:,3));

plot(tSol,ySol(:,4));

hold off

legend('x','vx','y','vy')

function dYdt=cometa(t,Y)

C=-1.33e20;

% x vx y vy % x=Y(1), vx=Y(2), y=Y(3), vy=Y(4)

r=sqrt((Y(1).^2+Y(3).^2));

dxdt=Y(2);

dvxdt=C.*Y(1)./r.^(3);

dydt=Y(4)

dvydt=C.*Y(3)./r.^(3);

dYdt=[dxdt;dvxdt;dydt;dvydt];

end

Answer 1

[Alan Stevens]

This works (I think your value of C is not in the right units):

tSpan=[0 2.33e9];
y0=[8.6e7 0 0 55.1];
[tSol,ySol]=ode45(@cometa,tSpan,y0);
plot(tSol,ySol(:,1));
hold on
plot(tSol,ySol(:,2));
plot(tSol,ySol(:,3));
plot(tSol,ySol(:,4));
legend('x','vx','y','vy')
hold off
figure(2)
plot(ySol(:,1),ySol(:,3),0,0,'*')
axis equal

function dYdt=cometa(~,Y)
M = 1.989*10^30; % kg
G = 6.674*10^-20; % km^3/(kg/s^2)
C=-M*G;
% x vx y vy  % x=Y(1), vx=Y(2), y=Y(3), vy=Y(4)
r=sqrt((Y(1).^2+Y(3).^2));
dxdt=Y(2);
dvxdt=C.*Y(1)./r.^3;
dydt=Y(4);
dvydt=C.*Y(3)./r.^3;
dYdt=[dxdt;dvxdt;dydt;dvydt];
end

[Stefano Rigante]

thank you very much for helping me.

i've tried your code , but when i run it and i plot only the value of tSol and ySol(:,1) i see an oscillating graph, but i expected a function with a max value.

[Alan Stevens]

You start the comet at x near perihelion with the sun to its "left". It will therefore move "left"; i.e. in a negative x direction (look at my figure(2)). The "maximum" x distance from the start will therefore be the negative of the minimum of x (you could also plot r as a function of time to see a maximum).. I suggest you plot the velocities on a separate graph from the distances as the scaling is very different.

[Stefano Rigante]

i've follwed your instruction , i've tried to find the minimum of x but it shows me -2.7258e5 that is minus than perihelion distance.

i've also tried to plot tSol and r as up defined but the graph that pops up is similar to the one oscillating of x.

can u show me how can i do that?

sorry for my incompetence

[Alan Stevens]

If you add

r = sqrt(ySol(:,1).^2 + ySol(:,3).^2);
figure(3)
plot(tSol,r),grid
xlabel('time'),ylabel('distance from sun')
disp(max(r))

just before the function definition you should see positive r, with a maximum.

Of course it cycles with time (your time span takes it just beyond one cycle).

[Stefano Rigante]

the max value of r i get is the initial value of x, that is 8.6e7km, same as periheliom . why??

[Alan Stevens]

I get 5.1619E9 km for aphelion. Run this code exactly as is to see:

tSpan=[0 2.33e9];
y0=[8.6e7 0 0 55.1];
[tSol,ySol]=ode45(@cometa,tSpan,y0);
plot(tSol,ySol(:,1));
hold on
plot(tSol,ySol(:,2));
plot(tSol,ySol(:,3));
plot(tSol,ySol(:,4));
legend('x','vx','y','vy')
hold off
figure(2)
plot(ySol(:,1),ySol(:,3),0,0,'*')
axis equal
r = sqrt(ySol(:,1).^2 + ySol(:,3).^2);
figure(3)
plot(tSol,r),grid
xlabel('time'),ylabel('distance from sun')
disp(max(r))
function dYdt=cometa(~,Y)
M = 1.989*10^30; % kg
G = 6.674*10^-20; % km^3/(kg/s^2)
C=-M*G;
% x vx y vy  % x=Y(1), vx=Y(2), y=Y(3), vy=Y(4)
r=sqrt((Y(1).^2+Y(3).^2));
dxdt=Y(2);
dvxdt=C.*Y(1)./r.^3;
dydt=Y(4);
dvydt=C.*Y(3)./r.^3;
dYdt=[dxdt;dvxdt;dydt;dvydt];
end

[Stefano Rigante]

You saved my life, thank u very much for the help.

the error was in the gravitational constant....

[Alan Stevens]

Numerical errors start to build up quickly as you do more and more cycles. That's probably why you seem to be getting decaying magnitudes (you might also get this with some of the other solvers).