MATLAB中求解分段函数定义的高度-雷诺数反转方程

问题背景 #

你需要反转「高度→雷诺数」的分段转换逻辑：已知马赫数（Ma）、特征长度（L）和雷诺数（Re），求解对应高度（h）。原转换函数以对流层顶阈值h_trop=11000m为分界，分两段计算Re，但直接用符号变量代入原函数求解时，会因h < h_trop这类条件判断无法作用于符号变量，出现「无法将符号变量转换为逻辑值」的错误。

解决方法 #

由于原函数是分段定义的，不能直接用符号变量代入整个函数求解，必须分区间单独建立符号方程，求解后验证结果是否落在对应区间内。具体步骤如下：

1. 明确两个区间的表达式 #

原函数分对流层（h < 11000m）和平流层下部（h ≥ 11000m）两段，我们需要分别写出两段的Re关于h的符号表达式，再解方程Re(h) = 目标Re。

2. 分区间求解并验证 #

• 先求解对流层区间的方程，检查结果是否小于11000m，符合则为有效解
• 再求解平流层区间的方程，检查结果是否大于等于11000m，符合则为有效解

完整MATLAB代码示例 #

% 已知参数
target_Re = 2.693844519000000e+07;
Ma = 0.78;
L = 4.371623535;

% 常量定义（和原函数一致）
h_trop = 11000;     % [m]
rho_0 = 1.225;      % [kg/m^3]
T_0 = 288.15;       % [K]
T_tropo = 216.65;   % [K]
T_s = 273.15;       % [K]
dT = -0.0065;       % [K/m]
k = 1.4;            % [-]
R = 287.058;        % [J/(K*kg)]
my_0 = 1.716e-005;  % [Pa*s]
C = 110.4;          % [K]
g = 9.80665;        % [m/s^2]
beta = my_0*(T_s+C)/(T_s^1.5);

% --------------------------
% 1. 求解对流层区间（h < 11000m）
% --------------------------
syms h_tropo real;
% 对流层的Re表达式
c_s_tropo = sqrt((T_0 + dT*h_tropo)*k*R);
v_tropo = Ma * c_s_tropo;
T_tropo_h = T_0 + dT*h_tropo;
my_tropo = beta*T_tropo_h^(1.5)/(T_tropo_h + C);
rho_tropo = rho_0*(1 + (dT*h_tropo/T_0))^(-g/R/dT - 1);
Re_tropo = rho_tropo * v_tropo * L / my_tropo;

% 建立方程并求解
eq_tropo = Re_tropo == target_Re;
sol_tropo = vpasolve(eq_tropo, h_tropo, [0, h_trop-1]);

% 验证解是否在区间内
if ~isempty(sol_tropo) && double(sol_tropo) < h_trop
    fprintf('对流层解：h = %.2f m\n', double(sol_tropo));
else
    fprintf('对流层区间无解\n');
end

% --------------------------
% 2. 求解平流层区间（h ≥ 11000m）
% --------------------------
syms h_strato real;
% 平流层的Re表达式
c_s_strato = sqrt(T_tropo*k*R);
v_strato = Ma * c_s_strato;
my_strato = beta*T_tropo^(1.5)/(T_tropo + C);
rho_strato = rho_0*(T_tropo/T_0)^(-g/R/dT - 1)*exp(-g/R/T_tropo*(h_strato - h_trop));
Re_strato = rho_strato * v_strato * L / my_strato;

% 建立方程并求解
eq_strato = Re_strato == target_Re;
sol_strato = vpasolve(eq_strato, h_strato, [h_trop, 20000]); % 平流层下部到20000m

% 验证解是否在区间内
if ~isempty(sol_strato) && double(sol_strato) >= h_trop
    fprintf('平流层解：h = %.2f m\n', double(sol_strato));
else
    fprintf('平流层区间无解\n');
end

补充说明 #

• 如果两个区间都得到有效解，需要结合实际飞行场景判断合理值（比如通常不会同时在对流层和平流层出现相同Re的情况）
• 如果两个区间都无解，说明输入的目标Re不在该Ma和L对应的高度范围内
• 求解时可以通过vpasolve的第三个参数指定搜索区间，缩小求解范围，提高效率和准确性

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