atan2 breaking ode15s

I have a program that runs ode15s a few thousand times in order to find a particular solution. However, I'm getting many integration tolerance errors such as the following:

"Warning: Failure at t=5.144337e+02.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.818989e-12) at time t."

Such warnings cause the program to slow down drastically, and sometimes even grind to a complete halt. The following is some test code that re-produces the error:

%Simulation constants
G = 6.672e-11; %Gravitational constant
M = 6.39e23; %Mass of Mars (kg)
g = 9.81; %Gravitational acceleration on Earth (m/s^2);
T1 = 845000/3; %Total engine thrust, 1 engine (N)
Isp = 282; %Engine specific impulse (s)
mdot1 = T1/(g*Isp); %Engine mass flow rate (kg/s)

xinit_on2 = [72368.8347685214;
            3384891.40103322;
            -598.36623436025;
            -1440.49702235844;
            16330.430678033]

tspan_on2 = [436.600093957202, 520.311296453027];

[t3,x3] = ode15s(@(t,x) engine_on_2(t, x, G, g, M, Isp, T1), tspan_on2, xinit_on2)

where the function engine_on_2 contains the system of ODEs that model the descent of a rocket, and is given by,

function xdot = engine_on_2(t, x, G, g, M, Isp, T1)
gamma = atan2(x(4),x(3)); %flight-path angle
xdot = [x(3); %xdot1: x-velocity
        x(4); %xdot2: y-velocity
        -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2))-(T1/x(5))*cos(gamma); %xdot3: x-acceleration
        -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2))-(T1/x(5))*sin(gamma); %xdot4: y-acceleration
        -T1/(g*Isp)]; %xdot5: engine mass flow rate
end

Having done some testing, it seems that I am getting the integration tolerance warnings because of the use of the atan2 function in

gamma = atan2(x(4),x(3))

which is used to calculate the flight-path angle of the rocket. If I change atan2 to another function (for example cos or sin) I don't get any integration tolerance warnings anymore (although, due to such a change, my solutions are obviously incorrect). As such, I was wondering if I am using atan2 incorrectly, or if there is a way to implement it differently so that I do not get the integration tolerance errors anymore. Furthermore, could it be that I am incorrect and that it is something other than atan2 that is causing the errors?

Answer 1

[Walter Roberson]

>> atan2(1e-5,1E-5)
ans =
         0.785398163397448
>> atan2(1e-5,-1E-5)
ans =
          2.35619449019234
>> atan2(-1e-5,1E-5)
ans =
        -0.785398163397448
>> atan2(-1e-5,-1E-5)
ans =
         -2.35619449019234

These possibilities represent trivial changes in position (on most scales), but get back quite different results. Therefore if you happen to cross one of the axes you could suddenly get big difference, which could give you problems.

For smoother calculations, you would probably be better replacing

gamma = atan2(x(4),x(3));
cos(gamma)
sin(gamma)

by

r = sqrt(x(3).^2 + x(4).^2);
cos_gamma = x(3) ./ r;
sin_gamma = x(4) ./ r;