MATLAB: atan2 breaking ode15s

Question

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?

Matt Matt · Accepted Answer · 2017-06-07T16:36:32

Use the odeset function to create an options structure that you then pass to the solver. RelTol and AbsTol can be adjusted in the ODE solver to eliminate your error. I was able to run your code using this addition without any errors:

options = odeset('RelTol',1e-13,'AbsTol',1e-20)

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

See the options are passed to the ODE solver as a 4th input parameter. Note the RelTol maxes out just above 1e-13 but hopefully that's fine for your application. Also you can try any of the other ODE solvers which can get rid of your error but from my playing around ode15s seems quite fast.

MATLAB: atan2 breaking ode15s

1 Answers