ode45 Matlab 中的 NaN 响应答案

【问题标题】：NaN response in ode45 Matlabode45 Matlab 中的 NaN 响应
【发布时间】：2018-09-19 00:46:22
【问题描述】：

我正在使用 ode45 来解决位置变量与时间的关系，但在初始条件之后，我只得到所有 4 个变量的 NaN 响应。这是我的功能代码：

function dxdt = ode(t,x)

global m b a x_teta k_lin k_tor Ip U Pinf gama 

m = 5;        % Mass
b = 1.5;      % Semi-chord
a = 0;        % Reference point to displacement -1 <= a <= 1
x_teta = 0;   % Static unbalanced parameter -> e - a
k_lin = 10;   % Linear stiffness
k_tor = 10;   % Torsional stiffness
Ip = 5;       % Inertia
U = 10;       % Velocity
Pinf = 1.184; % Air density

gama = 2*pi*Pinf*b*U^2; % Constant numbers from L 

% x(1) = h % x(2) = dh/dt % x(3) = teta % x(4) = dteta/dt

dxdt = [x(2);(-gama*x(3)-k_lin*x(1)-m*b*x_teta*(1/Ip)*(b*(.5+a)*gama*x(3)-k_tor*x(3)))/(m+(m^2*b^2*x_teta^2)/Ip);x(4);(-gama*x(3)-k_lin*x(1)-(b*(.5+a)*gama*x(3)-k_tor*x(3))/(b*x_teta))/(m*b*x_teta-(Ip)/(b*x_teta))]

以及实现代码：

tspan = 0:.01:20; x0 = [60; 0; 5; 0];

[t,x] = ode45(@ode,tspan,x0)

【问题讨论】：

标签： matlab ode45

【解决方案1】：

我怀疑设置这个 ODE 的物理有问题。

在 dxdt 的第 4 个元素的表达式中，有一个除以 x_teta，设置为零。

dxdt(4) = (-gama*x(3)-k_lin*x(1)-(b*(.5+a)*gama*x(3)-k_tor*x(3))/(b*x_teta))/(m*b*x_teta-(Ip)/(b*x_teta))

这个除以零导致一个NaN，并在迭代求解 ODE 时传播到其他元素。

【讨论】：