MATLAB：atan2 打破 ode15s答案

【问题标题】：MATLAB: atan2 breaking ode15sMATLAB：atan2 打破 ode15s
【发布时间】：2017-11-04 12:09:58
【问题描述】：

我有一个程序运行 ode15s 几千次以找到特定的解决方案。但是，我遇到了许多集成容差错误，例如：

"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."

此类警告会导致程序急剧变慢，有时甚至会完全停止。以下是一些重新产生错误的测试代码：

%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)

函数 engine_on_2 包含模拟火箭下降的 ODE 系统，由下式给出，

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

经过一些测试，我似乎收到了积分容差警告，因为在 gamma = atan2(x(4),x(3)) 中使用了 atan2 函数，该函数用于计算火箭的飞行路径角度。如果我将 atan2 更改为另一个函数（例如 cos 或 sin），我不会再收到任何积分容差警告（尽管由于这样的更改，我的解决方案显然是不正确的）。因此，我想知道我是否错误地使用了 atan2，或者是否有办法以不同的方式实现它，这样我就不会再遇到集成容差错误了。此外，可能是我不正确，并且导致错误的不是 atan2 吗？

【问题讨论】：

检查发布的解决方案，看看是否适合您。

标签： matlab numerical-integration

【解决方案1】：

使用odeset 函数创建一个选项结构，然后将其传递给求解器。 RelTol 和 AbsTol 可以在 ODE 求解器中进行调整以消除您的错误。我能够使用此添加运行您的代码而没有任何错误：

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)

请参阅options 作为第四个输入参数传递给 ODE 求解器。注意RelTol 的最大值刚好在 1e-13 以上，但希望这对您的应用程序来说没问题。您也可以尝试任何其他 ODE 求解器，它们可以消除您的错误，但从我玩 ode15s 来看似乎相当快。

【讨论】：