MATLAB 2022中使用ode45求解微分方程出现NaN结果问题

问题现象 #

使用ode45求解含流体阻力的运动微分方程时，除初始速度v0=0外，所有计算得到的v值均为NaN；仅将dvdt替换为2*t这类简单表达式时能得到有效结果。

核心问题分析 #

• 初始速度为0导致除以0错误
  当v=0时，雷诺数Re = (rho_W * v * D_B) / visc_W结果为0，计算阻力系数C_D时的24/Re会直接得到Inf，后续阻力F_D、加速度dvdt的计算都会因此变成NaN，导致ode45无法继续求解。
• 浮力与重力公式混淆
  原代码中浮力F_B的计算公式错误：浮力应为排开流体的重量，即F_B = rho_W * V_B * g，而非rho_B * V_B * g；而F_W = m_B * g = rho_B * V_B * g是物体自身重力，原代码中两者相等，导致F_B - F_W = 0，初始加速度仅由阻力项决定，进一步放大了NaN问题。

修复后的代码 #

clc
clear
% Define the time span
tspan = [0 10]; % Time span for the simulation (from 0 to 10 seconds)
v0 = 0;
% Solve the differential equation using ode45
[t, v] = ode45(@ode_function, tspan, v0);

% Plot the velocity as a function of time
plot(t, v)
xlabel('Time (s)');
ylabel('Velocity (m/s)');
title('Velocity vs. Time');
grid on;

function dvdt = ode_function(t, v)
    % Define parameters
    rho_B = 1000;     % Density of the body (kg/m^3)
    V_B = 0.01;       % Volume of the body (m^3)
    g = 9.81;         % Gravitational acceleration (m/s^2)
    rho_W = 1.225;    % Density of the fluid (kg/m^3)
    A_B = 0.1;        % Cross-sectional area of the body (m^2)
    m_B = rho_B * V_B; % Mass of the body (kg)
    D_B = 0.2;        % Diameter of the body (m)
    visc_W = 1.789e-5; % Dynamic viscosity of the fluid (N*s/m^2)
    
    % Handle v=0 case to avoid division by zero
    if v == 0
        % No drag force when velocity is 0
        F_D = 0;
        dvdt = (rho_W*V_B*g - m_B*g)/m_B;
        return;
    end
    
    % Calculate Reynolds number (Re)
    Re = (rho_W * v * D_B) / visc_W;
   
    % Calculate drag coefficient (C_D) using the given formula
    C_D = (24/Re) + (2.6 * (Re/5) / (1 + (Re/5)^1.52)) + (0.41 * (Re/263000)^-7.94 / (1 + (Re/263000)^-8)) + (Re^0.8 / 461000);
   
    % Calculate forces
    F_B = rho_W * V_B * g; % Correct buoyancy formula
    F_D = (0.5 * rho_W * A_B * C_D) * v^2;
    F_W = m_B * g;
   
    % Calculate acceleration (dv/dt)
    dvdt = (F_B - F_W - F_D)/m_B;
end

修复说明 #

1. 添加了v==0的判断分支，此时阻力F_D为0，直接计算无阻力时的加速度，避免除以0的错误。
2. 修正了浮力F_B的计算公式，确保浮力与重力的差值正确反映物体的受力趋势。

内容的提问来源于stack exchange，提问作者Ahmed Serag