MATLAB求解河流湍流扩散系数时数组C全为0的问题求助

问题描述 #

需要求解河流中的湍流扩散系数ε，因无法通过解析方法分离变量，故遍历ε∈[0,1]的数组，期望找到使浓度C(x,t)=3的ε值（t=1，其余参数为已知常数，另有两个变量对应数值表）。但运行代码后，C数组始终全为0。

原代码 #

clear
clc
Ma  =  9.35 * 10^12;
U = 93498.5;
x = 40233.6; 
t = 1;
epsilon =  linspace(0,1,100);
C = (Ma/(4*pi.*t.*epsilon)^0.5)*exp(-(x-(U.*t))^2/(4.*t.*epsilon));

问题原因分析 #

1. 指数项过小导致趋近于0：计算x - U*t得40233.6 - 93498.5 = -53264.9，平方后约为2.837×10^9。分母4*t*ε最大为4（当ε=1时），此时指数为-2.837×10^9 / 4 ≈ -7.09×10^8，exp(-7×10^8)的数值远小于MATLAB浮点数精度下限，直接返回0；当ε更小时，指数绝对值更大，结果同样为0。
2. ε=0的非法值：linspace(0,1,100)包含ε=0，此时分母(4*pi*t*epsilon)^0.5为0，虽理论上会得到无穷大，但结合后面的exp(-Inf)（0做分母导致指数项为负无穷），最终得到NaN，但因其他ε对应的C值全为0，整体数组显示为0。

解决方案 #

1. 检查参数合理性 #

确认x、U、t的单位是否统一，数值是否正确。当前U*t远大于x，意味着t=1时刻，污染物羽流中心还未到达x位置，浓度自然趋近于0。可尝试增大t值，使U*t接近x，缩小两者差值。

2. 修正代码避免非法值 #

将ε的起始值从0改为MATLAB最小正浮点数eps，避免除以0：

epsilon = linspace(eps, 1, 100);

3. 改用数值寻根方法 #

遍历数组精度有限且易错过解，推荐使用fzero函数直接寻找使C(ε)=3的解，同时先判断是否存在可行解：

clear
clc
Ma = 9.35e12;
U = 93498.5;
x = 40233.6; 
t = 1;

% 定义目标函数：C(epsilon) - 3 = 0
fun = @(epsilon) (Ma / sqrt(4*pi*t*epsilon)) * exp(-(x - U*t).^2/(4*t*epsilon)) - 3;

% 绘制C随ε的变化曲线，观察趋势
epsilon_vals = linspace(1e-6, 1, 100);
C_vals = fun(epsilon_vals) + 3;
plot(epsilon_vals, C_vals);
xlabel('湍流扩散系数ε');
ylabel('浓度C');
title('C随ε的变化趋势');

% 检查最大浓度是否达到3
max_C = max(C_vals);
disp(['当前参数下的最大浓度：', num2str(max_C)]);

if max_C < 3
    disp('当前参数组合下，不存在ε使C=3，请调整x、U、t或Ma的数值');
else
    % 寻找初始点并求解
    idx = find(C_vals >= 3, 1);
    epsilon_init = epsilon_vals(idx);
    epsilon_sol = fzero(fun, epsilon_init);
    disp('满足C=3的湍流扩散系数ε为：');
    disp(epsilon_sol);
end

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