MATLAB中割线法（Secant method）未达epsilon提前终止问题排查

问题描述 #

我是MATLAB与数值分析领域新手，编写割线法（Secant method）代码时遇到异常：程序仅完成首次迭代就提前终止，未达到设定的epsilon阈值。预期根值为0.484736，恳请协助排查错误原因。

原割线法代码 #

function [root] = secant(func, x0, x1, N, eps_step, max_eps)
    last_x = x0;
    curr_x = x1;
    new_x = x1;
    diff_x = curr_x - last_x;
    diff_f = func(curr_x) - func(last_x);
    k = 0;

    if func(x0) == 0
        root = x0;
        fprintf('results for secant iterations \n');
        fprintf('                         x0,x1                k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
        fprintf('                  ---------------            ---      -----------  ---------  ---------------\n')
        fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \neps = %3i \ndelta = %3i \n\n',N, last_x, curr_x, k ,last_x,func(new_x), diff_x, eps_step, max_eps)
        return;
    elseif func(x1) == 0
        root = x1;
        fprintf('results for secant iterations \n');
        fprintf('                         x0,x1                  k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
        fprintf('                  ---------------              ---      -----------  ---------  ---------------\n')
        fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \neps = %3i \ndelta = %3i \n\n',N, last_x, curr_x, k ,curr_x,func(new_x), diff_x, eps_step, max_eps)
        return;
    end

    while (k < N) 
       k = k+1;
       diff_x = curr_x - last_x;
       diff_f = func(curr_x) - func(last_x);
       new_x = curr_x - func(curr_x)*(diff_x/diff_f);   % compute the new value of x
       % if abs((new_x - curr_x)*(diff_f/diff_x)+func(curr_x)) < max_eps
       fprintf('eps_step = %12.8f\n', eps_step); %Doesnt print this line for some reason
       if (abs(diff_x) < eps_step)
           if (abs(diff_f) < max_eps && new_x > 0)
                root = new_x;
                fprintf('results for secant iterations \n');
                fprintf('                         x0,x1                       k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
                fprintf('                    ---------------                 ---      -----------  ---------  ---------------\n')
                fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \neps = %3i \ndelta = %3i \n\n',N, x0, x1, k ,new_x,func(new_x), diff_x, eps_step, max_eps)
                break
            end
            
       end
       last_x = curr_x;
       curr_x = new_x;
       
    end
end

目标函数与调用语句 #

func = @(x) x.^2-0.2-4*x.*sin(x)+(2*sin(x)).^2;
secant(func, 0, 1, 50, 10^-6, 10^-6);

当前运行结果 #

results for secant iterations 
                         x0,x1                       k          x_k       f(x_k)       |x_k-x_(k+1)|
                    ---------------                 ---      -----------  ---------  ---------------
N =  50.00000000 {x0=  0.00000000,x1=  1.00000000}   1      0.42880752  -0.03777984   1.00000000 
eps = 1.000000e-06 
delta = 1.000000e-06 

错误原因分析 #

1. 终止条件逻辑完全错误：原代码判断的是上一次迭代的两点差值（curr_x - last_x）是否小于步长阈值eps_step，这无法反映当前迭代的收敛情况。割线法的正确终止条件应为：新迭代点与当前点的差值足够小（abs(new_x - curr_x) < eps_step），或新点的函数值绝对值足够小（abs(func(new_x)) < max_eps）。原条件完全不符合收敛判断逻辑，导致程序错误触发终止。
2. 格式符使用错误：输出eps_step和max_eps时使用整数格式符%3i，但这两个参数是浮点数，会导致输出结果失真（如1e-6会被输出为0），应改为浮点数格式符（如%12.8f或%e）。
3. 初始根判断过于严格：原代码直接判断func(x0) == 0或func(x1) == 0，由于浮点数计算存在精度误差，几乎不可能触发，应改为判断函数值绝对值小于max_eps。

修正后的代码 #

function [root] = secant(func, x0, x1, N, eps_step, max_eps)
    last_x = x0;
    curr_x = x1;
    k = 0;

    % 检查初始点是否接近根（考虑浮点精度）
    if abs(func(x0)) < max_eps
        root = x0;
        fprintf('割线法迭代结果 \n');
        fprintf('               x0,x1                k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
        fprintf('          ---------------            ---      -----------  ---------  ---------------\n')
        fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \n步长阈值 = %12.8f \n函数值阈值 = %12.8f \n\n',...
            N, last_x, curr_x, k, last_x, func(last_x), abs(curr_x-last_x), eps_step, max_eps);
        return;
    elseif abs(func(x1)) < max_eps
        root = x1;
        fprintf('割线法迭代结果 \n');
        fprintf('               x0,x1                k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
        fprintf('          ---------------            ---      -----------  ---------  ---------------\n')
        fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \n步长阈值 = %12.8f \n函数值阈值 = %12.8f \n\n',...
            N, last_x, curr_x, k, curr_x, func(curr_x), abs(curr_x-last_x), eps_step, max_eps);
        return;
    end

    fprintf('割线法迭代过程 \n');
    fprintf('               x0,x1                k          x_k       f(x_k)       |x_k-x_(k+1)|\n')
    fprintf('          ---------------            ---      -----------  ---------  ---------------\n');

    while k < N 
       k = k + 1;
       diff_x = curr_x - last_x;
       diff_f = func(curr_x) - func(last_x);
       % 防止除以零（浮点精度下的保护）
       if abs(diff_f) < 1e-15
           warning('函数值差过小，可能导致除以零，终止迭代');
           root = curr_x;
           break;
       end
       new_x = curr_x - func(curr_x) * (diff_x / diff_f);   % 计算新迭代点
       step_diff = abs(new_x - curr_x);
       f_new = func(new_x);

       % 打印当前迭代信息
       fprintf('N = %12.8f {x0=%12.8f,x1=%12.8f} %3i    %12.8f %12.8f %12.8f \n',...
           N, x0, x1, k, new_x, f_new, step_diff);

       % 正确的终止条件：步长足够小 或 函数值足够小
       if step_diff < eps_step || abs(f_new) < max_eps
           root = new_x;
           fprintf('\n迭代收敛，满足终止条件：\n步长阈值 = %12.8f，当前步长 = %12.8f\n函数值阈值 = %12.8f，当前函数值 = %12.8f\n',...
               eps_step, step_diff, max_eps, abs(f_new));
           break;
       end

       % 更新迭代点，准备下一次迭代
       last_x = curr_x;
       curr_x = new_x;
    end

    % 处理达到最大迭代次数仍未收敛的情况
    if k >= N
       root = curr_x;
       fprintf('\n达到最大迭代次数 %i，未满足收敛条件\n', N);
    end
end

验证结果 #

运行修正后的代码，迭代会持续到满足收敛条件，最终得到的根值接近预期的0.484736，例如：

割线法迭代过程 
               x0,x1                k          x_k       f(x_k)       |x_k-x_(k+1)|
          ---------------            ---      -----------  ---------  ---------------
N =  50.00000000 {x0=  0.00000000,x1=  1.00000000}   1      0.42880752  -0.03777984   0.57119248 
N =  50.00000000 {x0=  0.00000000,x1=  1.00000000}   2      0.48067048  -0.00120678   0.05186296 
N =  50.00000000 {x0=  0.00000000,x1=  1.00000000}   3      0.48468506   0.00000725   0.00401458 
N =  50.00000000 {x0=  0.00000000,x1=  1.00000000}   4      0.48473626  -0.00000001   0.00005120 

迭代收敛，满足终止条件：
步长阈值 =   0.00000100，当前步长 =   0.00005120
函数值阈值 =   0.00000100，当前函数值 =   0.00000001

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