在MATLAB 2022中使用ode45时的未定义结果(NaN)

cwtwac6a  于 2023-10-23  发布在  Matlab
关注(0)|答案(1)|浏览(351)

验证码:

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)
    % 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_B * V_B * g;
    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_D - F_W))/m_B);
end

这段代码应该在使用时输出一个正确的图,但是v的所有值(除了第一个值为0)都是undefined/NaN。这个问题可能是在ode_function()函数中,但我对ode 45没有太多经验。
我试着修改dvdt的定义,当把它变成一个更简单的方程,比如2t时,得到了正确的结果,但这是我能得到的有用的结果。

hvvq6cgz

hvvq6cgz1#

通过使用断点进行调试(您应该自己检查debugging in Matlab),我发现了以下事情:
初始化v0=0并在ode_function中将其用作Re = (rho_W * v * D_B) / visc_W的因子,因此是Re=0
这将导致公式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),您将其除以(24/Re-> Inf)。在公式的后面,你可以提高Re(其中一种情况是Re^0.8-> Inf)。这导致Inf/Inf-> NaN。
由于这一点,你做的任何事情都不会产生任何有用的结果,包括情节。
简单的解决方案是不使用0作为输入v0和/或验证(作为<>0)它作为函数的参数

相关问题