基于求积法的while无限循环问题调试求助

问题概述 #

这段Matlab/Octave代码存在两个核心问题：

• 使用while n >= 1作为循环条件时，程序陷入无限循环，因为终止循环的if条件始终无法满足
• 若将循环条件改为有限范围（比如while n <= 100），循环能正常终止，但if语句仍不会执行，最终导致n_escolhido变量未定义

核心问题分析 #

1. 导数计算完全错误
   原函数y(t) = t*exp(-t/4)的导数需用乘积法则计算：
   y’(t) = exp(-t/4) + t*(-1/4)*exp(-t/4) = exp(-t/4)*(1 - t/4)
   但代码中y_deriv的定义写成了((-(t-4)/4) * e^(t/4))^2，不仅符号、指数项错误，还额外添加了平方，直接导致F(z)的计算值完全偏离正确结果，后续积分差值永远达不到终止阈值。

2. Simpson积分实现逻辑混乱

• 内层嵌套的while k <=n-1和for l=1:length(T)完全多余，T是长度为31的向量，这会让每个积分区间被重复计算31次，导致求和结果严重偏大，积分值完全错误。
• Qn和Q2n的计算方式错误：Simpson积分要求Qn是n段（偶数段）的积分值，Q2n是2n段的积分值，代码中只是简单将系数除以2，并未重新计算2n段的积分。
• valor_3错误赋值为F(xi_2)，应该对应F(xi_3)。

3. 阈值写法错误
   代码中10e-6实际等于1e-5，如果要表示10^-6，应该写成1e-6。

修复后的代码 #

% 定义参数与函数
x = @(t) t;
y = @(t) t*exp(-t/4);
a = 0;
b = 15;

% 定义导数（修正错误）
x_deriv = 1;
y_deriv = @(t) exp(-t/4) * (1 - t/4);
% 定义弧长积分函数
F = @(z) sqrt(x_deriv^2 + y_deriv(z)^2);

n = 1;
max_n = 100; % 设置循环上限，防止极端情况无限循环
tolerance = 1e-6; % 修正阈值为10^-6
n_escolhido = [];

while n <= max_n
    % 计算Qn：n段Simpson积分（n需为偶数，这里自动调整为最近的偶数）
    n_even = 2 * ceil(n/2);
    h_n = (b - a)/n_even;
    t_n = a:h_n:b;
    f_n = F(t_n);
    Qn = h_n/3 * (f_n(1) + f_n(end) + 4*sum(f_n(2:2:end-1)) + 2*sum(f_n(3:2:end-2)));
    
    % 计算Q2n：2n段Simpson积分
    n2_even = 2 * n_even;
    h_2n = (b - a)/n2_even;
    t_2n = a:h_2n:b;
    f_2n = F(t_2n);
    Q2n = h_2n/3 * (f_2n(1) + f_2n(end) + 4*sum(f_2n(2:2:end-1)) + 2*sum(f_2n(3:2:end-2)));
    
    % 检查终止条件
    if abs(Qn - Q2n) <= tolerance
        n_escolhido = n_even;
        break;
    end
    
    n = n + 1;
end

if ~isempty(n_escolhido)
    disp(['满足精度要求的n值：', num2str(n_escolhido)]);
else
    disp(['达到最大循环次数', num2str(max_n), '仍未满足精度要求']);
end

关键修改说明 #

• 修正导数计算：还原y_deriv的正确表达式，确保F(z)计算的是弧长的被积函数。
• 重构Simpson积分逻辑：分别实现n段和2n段的Simpson积分计算，遵循标准的Simpson公式：∫f(x)dx ≈ h/3*(f0 + 4f1 + 2f2 + 4f3 + ... + 4f(n-1) + fn)。
• 调整循环逻辑：设置循环上限max_n，避免极端情况下的无限循环；自动将n调整为偶数（Simpson积分要求分段数为偶数）。
• 修正阈值：将10e-6改为1e-6，匹配需求中的10^-6精度。
• 增加异常处理：判断n_escolhido是否为空，避免未定义变量报错。

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