Thread Subject: ode45 changes solution with units of measuremetns

Hi everyone!

I've got a strange problem with ode45. My task is rather simple: there is equations of motion of particle in the galactic gravitational potential, which must be integrated over timespan of 4.5 Gyr. When I change units of measurements for any parameters (for example, time), ode45 gives solution very different from the previous. Here is the function, motion10, with comments on changes in function motion8. Motion10 and motion8 uses different units of time (8e4 and 5e4 yrs), which gives different solution for equal initial conditions. Of course, all constants have been changed according to change of time units.

function HVeloc = motion10(Nh)

% UNITS:
% [R] kpc, [t] 8e4 yrs in motion10, 5e4 yrs in motion8

VtoRt = 8.1761e-5; % [km/s] to [kpc/8e4 yrs] {[kpc/5e4 yrs]} in motion10 {motion8}
RttoV = 1/VtoRt; % VtoRt = 5.1101e-5 in motion8

% Gravitational constant in [Msun,kpc,8e4 yrs] {[Msun,kpc,5e4 yrs]}
G = 1.5625*4.94e-14; % or 4*4.94e-14 in motion8

% Potential param-s in [Msun, kpc]
a = 1e-3*[0,3.7e3];
b = 1e-3*[0.277e3,0.2e3];
M = [1.12e10,8.07e10,5e10];
rh = 1e-3*6e3;

% initial position and velocity
R0 = 100; % in kpc
z0 = 100;
    
VR0 = 20; % in km/s
Vz0 = 20;
Vt0 = 20;
    
tspan=[0 56250]; % tspan (4.5e9 yrs) in units of 8e4 yrs, [0 90000] in motion8

% velocity from [km/s] to [kpc/8e4 yrs]
R0dyn = R0;
z0dyn = z0;
VR0dyn= VtoRt*VR0;
Vz0dyn= VtoRt*Vz0;
Vt0dyn= VtoRt*Vt0;

L = Vt0dyn*R0dyn;
    
% solving ODE
options = odeset('NonNegative',1,'RelTol',1e-10);
y0=[R0dyn; z0dyn; VR0dyn; Vz0dyn];
[t,y] = ode45(@motode, tspan, y0,options); % y=[Rdyn;zdyn;VRdyn;Vzdyn]
    
TT = size(t);
    
figure(1);
plot(y(:,1),y(:,2));
hold on;
    
Vh = sqrt(y(TT(1),3).^2 + y(TT(1),4).^2); % velocity at the end
Vn = RttoV*Vh; % in units of [km/s]
HVeloc = Vn; % result

% potential and its derivatives
function U = Ugrav (R,z)
    
    U = -G*M(1)/(R^2+(a(1)+(z^2+b(1)^2)^(1/2))^2)^(1/2)-G*M(2)/(R^2+(a(2)+(z^+b(2)^2)^(1/2))^2)^(1/2)-G*M(3)/rh*(1/2*log(1+(R^2+z^2)/rh^2)+rh/(R^2+z^2)^(1/2)*atan((R^2+z^2)^(1/2)/rh));

end

function dU = dUR (R,z)

    dU = G*M(1)/(R^2+(a(1)+(z^2+b(1)^2)^(1/2))^2)^(3/2)*R+G*M(2)/(R^2+(a(2)+(z^2+b(2)^2)^(1/2))^2)^(3/2)*R-G*M(3)/rh*(R/rh^2/(1+(R^2+z^2)/rh^2)-rh/(R^2+z^2)^(3/2)*atan((R^2+z^2)^(1/2)/rh)*R+1/(R^2+z^2)*R/(1+(R^2+z^2)/rh^2));

end

function dU = dUz (R,z)

    dU = G*M(1)/(R^2+(a(1)+(z^2+b(1)^2)^(1/2))^2)^(3/2)*(a(1)+(z^2+b(1)^2)^(1/2))/(z^2+b(1)^2)^(1/2)*z+G*M(2)/(R^2+(a(2)+(z^2+b(2)^2)^(1/2))^2)^(3/2)*(a(2)+(z^2+b(2)^2)^(1/2))/(z^2+b(2)^2)^(1/2)*z-G*M(3)/rh*(z/rh^2/(1+(R^2+z^2)/rh^2)-rh/(R^2+z^2)^(3/2)*atan((R^2+z^2)^(1/2)/rh)*z+1/(R^2+z^2)*z/(1+(R^2+z^2)/rh^2));

end

% pass to ode45
    function dRzdt = motode(t,y)
        
        dRzdt = [ y(3);
                           y(4);
                dUR(y(1),y(2)) + L^2/y(1)^3;
                dUz(y(1),y(2)) ];
    end

end